W badaniu wykorzystano publicznie dostępne zestawy danych oraz komercyjnie dostępne ludzkie linie komórkowe. Nie brali w nim udziału uczestnicy będący ludźmi, nie wykorzystano nowo pobranych tkanek ludzkich ani kręgowców; w związku z tym nie wymagano zgody instytucjonalnej komisji etyki. Narzędzia badawcze, oprogramowanie i zasoby online wykorzystane w protokole wymieniono w Tabeli materiałów.
1. Przetwarzanie danych dla analizy ekspresji różnicowej
W celu stworzenia kohort odkrycia i walidacji transkrypcyjnej przeszukano bazę danych Gene Expression Omnibus (GEO) przy użyciu terminu „ulcerative colitis”. Cztery zbiory danych — GSE87466, GSE92415, GSE107499 oraz GSE75214 — posłużyły jako zestaw treningowy, natomiast GSE47908 zarezerwowano do walidacji zewnętrznej. Zbiory danych treningowych znormalizowano i poddano korekcie efektu serii za pomocą narzędzia ComBat w programie Sangerbox 3.0, natomiast zbiór walidacyjny przetworzono niezależnie.
Łącznie 583 genów związanych z ferroptozą (FRGs) pozyskano z bazy FerrDB, a następnie uzupełniono o geny opisane w literaturze w celu zredukowania błędu wynikającego z korzystania z jednego źródła12. Zbiór 96 genów związanych z kuproptozą (CRGs) zestawiono na podstawie opublikowanych badań12,13,14. Analizę różnicowej ekspresji przeprowadzono przy użyciu pakietu limma, przyjmując progi |log₂ fold change| > 1 oraz wartość P skorygowaną metodą Benjamini-Hochberga < 0,05. Wykresy diagnostyczne normalizacji, korekcji efektu serii, różnicowej ekspresji oraz klasteryzacji przedstawiono na Rysunku uzupełniającym 1. Skład próbek, wykorzystane platformy oraz przydziały kohort podano w Tabeli uzupełniającej 1.
2. Identyfikacja CF-DEG w UC
Analizę korelacji Pearsona przeprowadzono przy zastosowaniu progów |r| > 0,5 i P < 0,05 w celu zidentyfikowania współwyrażonych genów związanych z ferroptozą (FRGs) oraz genów związanych z kuproptozą (CRGs). Otrzymany zestaw genów zestawiono z genami o różnej ekspresji (DEGs), aby zidentyfikować współwyrażone geny o różnej ekspresji związane z kuproptozą i ferroptozą (CF-DEGs).
3. Analiza wzbogaceń
Analizy wzbogacenia według Gene Ontology (GO) oraz Kyoto Encyclopedia of Genes and Genomes (KEGG) przeprowadzono przy użyciu narzędzia clusterProfiler oraz platformy Bioinformatics.com.cn. Jako próg istotności przyjęto wartość P < 0,05 po korekcie metodą Benjamini–Hochberga. Wykresy wzbogacenia przedstawiono na Rysunku uzupełniającym 2.
4. Ważona analiza sieci koekspresji genów (WGCNA)
Strukturę modularną leżącą u podstaw profili ekspresji genów przeanalizowano za pomocą ważonej analizy sieci współekspresji genów (WGCNA). Próbki nieprawidłowe zidentyfikowano i wykluczono poprzez grupowanie hierarchiczne w połączeniu z funkcją goodSamplesGenes. Następnie wyznaczono odpowiednią wartość potęgi miękkiego progowania (β), aby spełnić kryterium topologii wolnej od skali niezbędne do konstrukcji sieci.
Przetworzone dane dotyczące ekspresji wykorzystano do skonstruowania macierzy nakładania topologicznego (TOM) na potrzeby analizy koekspresji genów. Z zastosowaniem metody dynamicznego cięcia drzewa geny o wysoce spójnych wzorcach ekspresji zorganizowano w moduły genowe, z których każdy zawierał co najmniej 100 genów. Powiązania między eigengenami modułów (ME) a fenotypami przeanalizowano w celu zidentyfikowania modułów związanych z procesem chorobowym i wyodrębnienia odpowiadających im zestawów genów. Diagnostyka przesiewowa WGCNA została przedstawiona na Rysunku uzupełniającym 3.
5. Wybór biomarkerów
Aby wyłonić priorytetowe kandydaty do dalszych badań eksperymentalnych, równolegle zastosowano trzy algorytmy uczenia maszynowego15. Regresja LASSO (least absolute shrinkage and selection operator) została wykorzystana do selekcji cech, las losowy (random forest) z ntree = 500 posłużył do określenia rangi istotności genów, a maszyna wektorów wspierających z rekurencyjną eliminacją cech (SVM-RFE) została użyta do zidentyfikowania podzbioru cech z najniższym błędem walidacji krzyżowej. Pełna diagnostyka przesiewowa znajduje się na Rysunku uzupełniającym 5.
Znaczenie genów oceniano na podstawie wkładu poszczególnych genów w dokładność klasyfikacji, a geny o najwyższej randze zostały wyselekcjonowane w celu stworzenia wstępnego zestawu genów16. Algorytm SVM-RFE przeprowadził przesiewanie cech i optymalizację modelu poprzez iteracyjną eliminację cech, co pozwoliło na identyfikację optymalnego podzbioru cech dla klasyfikacji17. Geny wspólne dla trzech algorytmów uznano za potencjalne biomarkery zidentyfikowane za pomocą uczenia maszynowego.
Wykorzystano CF-DEG do skonstruowania sieci oddziaływań białko-białko (PPI), a następnie użyto narzędzia CytoHubba do rankingowania węzłów kluczowych. Przecięcie kluczowych modułów WGCNA, trzech wyników uczenia maszynowego oraz zestawu rdzeniowego PPI pozwoliło zdefiniować pięć zintegrowanych potencjalnych biomarkerów: LCN2, IDO1, CXCL2, NOS2 oraz CD274. Geny te sklasyfikowano jako diagnostyczne kandydatów związanych z UC. W dalszych eksperymentach z inhibitorami LCN2, IDO1, CXCL2 i NOS2 sklasyfikowano jako markery reagujące na hamowanie ferroptozy, natomiast CD274 jako marker reagujący na hamowanie kuproptozy. Geny te nie zostały uznane za udowodnionych regulatorów nadrzędnych żadnej z tych ścieżek.
6. Opracowanie modelu prognostycznego i ewaluacja wielometodowa
Aby ocenić potencjał translacyjny w zastosowaniach klinicznych zidentyfikowanej sygnatury genów, przy użyciu pakietu rms skonstruowano nomogram służący do oceny skuteczności diagnostycznej. Wydajność modelu oceniono za pomocą krzywej kalibracyjnej, testu Hosmera–Lemeshow (HL) oraz indeksu zgodności Harrella (C-index).
Kalibrację modelu oceniono za pomocą testu HL oraz średniego błędu bezwzględnego (MAE). Przyjęto, że wartość P > 0,05 oraz MAE < 0,1 będą odpowiednio progami akceptowalnego dopasowania i dokładności modelu. Wygenerowano krzywe charakterystyki operacyjnej odbiornika (ROC) oraz obliczono pole pod krzywą (AUC), aby ocenić wydajność prognostyczną nomogramu dla UC. Do oceny generalizacji i stabilności modelu wykorzystano niezależny zbiór danych.
7. Analiza wzbogacenia zbiorów genów (GSEA)
W zbiorze treningowym przeprowadzono analizę korelacji Spearmana pomiędzy każdym biomarkerem a wszystkimi pozostałymi genami. Wyniki korelacji zostały uszeregowane w celu wygenerowania uporządkowanego zbioru genów. Analizę wzbogacenia zbiorów genów (GSEA) przeprowadzono z wykorzystaniem zbioru genów c2.cp.kegg_medicus.v2025.1.Hs.symbols z bazy Molecular Signatures Database (MSigDB) jako referencyjnego. Ścieżki biologiczne z P < 0,05 uznano za istotnie wzbogacone.
8. Analiza krajobrazu immunologicznego
Narzędzie CIBERSORT oraz macierz sygnatur LM22 zostały wykorzystane do przeprowadzenia analizy dekonwolucji danych ekspresyjnych i oszacowania proporcji 22 podtypów komórek odpornościowych. Następnie przeprowadzono analizę korelacji Spearmana w celu oceny powiązań między pięcioma biomarkerami a podtypami komórek odpornościowych.
9. Analiza sieci regulacyjnych
Informacje na temat regulacji przez czynniki transkrypcyjne (TF) oraz mikroRNA (miRNA) powiązane z genami uzyskano z zintegrowanych baz danych ChEA3 i TarBase 9.0 za pośrednictwem platformy NetworkAnalyst. Do sieci regulacyjnej włączono miRNA powiązane z co najmniej dwoma węzłami.
10. Sekwencjonowanie RNA z pojedynczych komórek (Single-cell RNA-seq)
Zintegrowano zbiory danych z pojedynczych komórek dla UC: GSE116222, obejmujący trzy próbki zdrowych i trzy próbki z UC, oraz GSE214695, obejmujący sześć próbek zdrowych i sześć próbek z UC. Zachowano komórki, w których wykryto 200–5 000 genów, >1 000 unikalnych identyfikatorów molekularnych (UMI), < 30% transkryptów mitochondrialnych, < 3% transkryptów hemoglobiny oraz < 50% transkryptów rybosomalnych. Do normalizacji wykorzystano funkcję NormalizeData, a do korekty efektu seryjnego przed klastrowaniem i adnotacją zastosowano narzędzie Harmony.
Po filtrowaniu dane zostały znormalizowane przy użyciu funkcji NormalizeData. Do minimalizacji efektu serii pomiędzy próbkami wykorzystano algorytm Harmony. W celu redukcji wymiarowości i klastrowania komórek, aby scharakteryzować struktury lokalne i globalne w obrębie populacji komórek, zastosowano metodę UMAP (Uniform Manifold Approximation and Projection) oraz t-SNE (t-distributed Stochastic Neighbor Embedding). Do dalszych analiz zachowano łącznie 31 712 wysokiej jakości komórek. Typy komórek przypisano na podstawie kanonicznych genów markerowych dla poszczególnych typów komórek.
Analizę czasu pseudoczynnościowego (pseudotime) przeprowadzono przy użyciu oprogramowania Monocle 2 w celu zbadania zmian w stanach komórek nabłonkowych i mieloidalnych z tkanek UC. Ewaluowano trajektorie rozwoju komórek oraz różnice między grupami, a następnie zwizualizowano wzorce ekspresji biomarkerów wzdłuż trajektorii czasu pseudoczynnościowego.
Do oceny komunikacji międzykomórkowej pomiędzy subpopulacjami komórek nabłonkowych i mieloidalnych wykorzystano narzędzie CellChat. Zintegrowano wzorce ekspresji specyficzne dla poszczególnych subpopulacji, aby wywnioskować prawdopodobieństwa oddziaływań oraz kluczowe szlaki sygnalizacyjne między komórkami nabłonkowymi a mieloidalnymi i scharakteryzować ich potencjalne sieci regulacji międzykomórkowej.
11. Eksperymenty komórkowe
Ludzkie komórki raka gruczołowego jelita grubego Caco-2 hodowano w pożywce Dulbecco’s modified Eagle zawierającej 10% płodowej surowicy bydlęcej, 100 U/mL penicyliny i 100 µg/mL streptomycyny w temperaturze 37 °C w atmosferze 5% CO₂. Komórki Caco-2 wybrano jako powtarzalny model nabłonkowej monowarstwy do testów bariery jelitowej i uszkodzeń zapalnych; jednak ich nowotworowe pochodzenie ogranicza bezpośrednią generalizację wyników na nienowotworowy nabłonek jelita grubego. Informacje o pasażowaniu nie były dostępne w archiwalnych zapisach eksperymentalnych.
Model stanu zapalnego i uszkodzenia in vitro opracowano poprzez narażenie komórek na 1 µg/mL lipopolisarydu (LPS) przez 24 h. Komórki przydzielono do grupy kontrolnej z nośnikiem w postaci dimetylosulfoksydu, grupy z samym LPS, grupy z LPS + 4 µM ferrostatiny-1 lub grupy z LPS + 40 µM tetratiomolibdanu, przy zachowaniu dopasowanych końcowych stężeń rozpuszczalnika. Grupa z samym LPS służyła jako pozytywna kontrola uszkodzenia/modelu do porównań z grupami z zastosowaniem inhibitorów.
- Żywotność komórek
Żywotność komórek oceniono za pomocą testu żywotności komórek zgodnie z instrukcjami producenta. Po inkubacji w temperaturze 37 °C przez około 2 h zmierzono absorbancję przy 450 nm za pomocą czytnika mikropłytek.
- Pomiar poziomów MDA i Fe²⁺
Wewnątrzkomórkowe poziomy malondialdehydu (MDA) i Fe2⁺ zmierzono przy użyciu odpowiednich zestawów testowych zgodnie z instrukcjami producentów. Całkowite białko oznaczono metodą kwasu bicinchoninowego, a każdy odczyt znormalizowano do zawartości białka w odpowiedniej próbce.
- Analiza z użyciem fluorescencyjnej sondy jonów miedzi
Po usunięciu pożywki hodowlanej dodano specyficzną dla Cu2⁺ sondę fluorescencyjną w końcowym stężeniu 1 µM i inkubowano przez 30 min w temperaturze 37 °C. Obrazy pozyskano natychmiast, bez przemywania, za pomocą odwróconego mikroskopu fluorescencyjnego przy długościach fal wzbudzenia i emisji odpowiednio 510 i 578 nm. Parametry akwizycji pozostały stałe dla wszystkich grup, a skala obrazu została skalibrowana na 50 µm.
- RT-qPCR
Całkowite RNA wyizolowano przy użyciu odczynnika do ekstrakcji RNA opartego na fenolu i guanidynium. Komplementarne DNA zsyntetyzowano z 1 µg RNA przy użyciu zestawu do syntezy pierwszej nici cDNA, natomiast do reakcji PCR w czasie rzeczywistym wykorzystano odczynniki oparte na SYBR Green. Każdą reakcję przeprowadzono w trzech powtórzeniach technicznych. Sekwencje starterów podano w Tabeli uzupełniającej 2, a względną obfitośćt mRNA obliczono metodą 2⁻ΔΔCt.
12. Analiza statystyczna
Normalność rozkładu oraz jednorodność wariancji oceniono odpowiednio za pomocą testów Shapiro–Wilka oraz Levene'a z centrowaniem do mediany. Dane spełniające oba założenia analizowano za pomocą nieparzystego dwustronnego testu t lub jednoczynnikowej analizy wariancji, a następnie testu Tukey’a (HSD). W przypadku niespełnienia któregokolwiek z założeń stosowano test t Welcha lub analizę wariancji Welcha, a następnie porównania Gamesa–Howella.
Wyniki przedstawiono jako średnia ± odchylenie standardowe. Rzeczywista liczba powtórzeń biologicznych, dokładne dwustronne wartości P oraz weryfikacje założeń zostały zamieszczone w Tabeli uzupełniającej 3. Zarchiwizowane arkusze źródłowe dla Ryciny 13 zawierają mniejszą liczbę powtórzeń kontrolnych dla paneli A (n = 4), C (n = 4) i F (n = 2), podczas gdy wszystkie pozostałe kombinacje paneli/grup obejmują n = 6.
Analizy obliczeniowe przeprowadzono przy użyciu programu R oraz narzędzi limma, WGCNA, clusterProfiler, glmnet, randomForest, e1071, pROC, rms, Seurat v4, Harmony, Monocle 2, CellChat oraz Cytoscape/CytoHubba. Wykorzystano zasoby internetowe, takie jak FerrDB, GEO, STRING, ChEA3, TarBase 9.0, NetworkAnalyst, CIBERSORT/LM22, MSigDB, Bioinformatics.com.cn oraz Sangerbox 3.0.