Oświadczenie komisji etyki
Badanie to zostało przeprowadzone zgodnie z Deklaracją Helsińską. Protokół został zatwierdzony przez Komisję Etyczną Szpitala Tradycyjnej Medycyny Chińskiej Shenzhen Luohu (numer zatwierdzenia 2024-LHQZYYYXLL-KY-039), a od wszystkich uczestników przed włączeniem do badania uzyskano pisemną świadomą zgodę. Szczegółowe informacje na temat narzędzi badawczych i materiałów wykorzystanych w tym protokole znajdują się w Tabeli materiałów.
Źródło i przetwarzanie danych
Zbiory danych dotyczących ekspresji genów powiązanych z POChP pobrano z bazy Gene Expression Omnibus (GEO). Zbiór GSE54837 wykorzystano jako zbiór transkryptomiczny, a zbiór GSE112811 posłużył jako zestaw walidacyjny (Tabela 1). Geny ac4C-RG zgromadzono na podstawie literatury18. Geny różnicowo ulegające ekspresji (DEG) pomiędzy grupą z POChP a grupą kontrolną zidentyfikowano przy użyciu pakietu R limma. Geny DEG uznano za statystycznie istotne, jeśli |log2FC| > 0 oraz p < 0,05. W celu wizualizacji ogólnego rozkładu zmian ekspresji genów wygenerowano wykresy wulkaniczne.
Konstrukcja WGCNA
Analizę WGCNA przeprowadzono na zbiorze danych GSE54837 przy użyciu języka R w celu zidentyfikowania modułów związanych z POChS. Przed konstrukcją sieci zidentyfikowano i usunięto próbki odstające za pomocą hierarchicznej analizy skupień, wykorzystując funkcję hclust z metodą średniego wiązania (average linkage) oraz miarą odległości euklidesowej. Optymalna wartość parametru miękkiego progowania (soft-thresholding power) (β = 10) wybrano, aby uzyskać wskaźnik dopasowania topologii bezskalowej R2 ≥ 0,85, zapewniając równowagę między topologią wolną od skali a średnią łącznością. Skonstruowano macierz sąsiedztwa, którą następnie przekształcono w macierz nakładania topologicznego (TOM). Moduły genów zidentyfikowano przy użyciu algorytmu dynamicznego cięcia drzewa (deepSplit = 2, minClusterSize = 50). Moduły o korelacjach eigen-genów > Następnie, przy użyciu funkcji mergeCloseModules, scalono moduły o wartości 0,75. W celu zidentyfikowania modułów związanych z POChP do dalszej analizy, eigengeny modułów skorelowano z cechami klinicznymi (status POChP, wiek, płeć i status palenia) przy użyciu współczynników korelacji Pearsona.
Przesiewanie, analiza wzbogacenia i analiza sieci PPI dla genów nakładających się
W celu zidentyfikowania genów wspólnych dla DEG, genów modułu MEsalmon oraz ac4C-RG wygenerowano wykres Venna przy użyciu pakietu R ggvenn. Analizę wzbogacenia funkcjonalnego genów wspólnych przeprowadzono z wykorzystaniem baz danych Gene Ontology (GO) oraz Kyoto Encyclopedia of Genes and Genomes (KEGG) przy użyciu pakietu R clusterProfiler. Informacje o oddziaływaniach białko-białko (PPI) pobrano z bazy danych STRING (https://string-db.org/), aby przeanalizować interakcje na poziomie białek pomiędzy genami wspólnymi. Do wizualizacji wynikowej sieci PPI wykorzystano oprogramowanie Cytoscape.
Identyfikacja kluczowych genów za pomocą uczenia maszynowego
Zastosowano trzy techniki uczenia maszynowego: regresję LASSO (least absolute shrinkage and selection operator), XGBoost (extreme gradient boosting) oraz las losowy (RF). Regresję LASSO zaimplementowano przy użyciu pakietu glmnet z 10-krotną walidacją krzyżową w celu wyznaczenia optymalnego parametru kary λ. Parametr type.measure ustawiono na "deviance", a parametr family na "binomial". Optymalne λ wybrano na podstawie kryterium λmin, które minimalizuje dewiancję z walidacji krzyżowej, co pozwoliło wyłonić 17 genów. XGBoost wykonano przy użyciu pakietu xgboost z następującymi hiperparametrami: nrounds = 100, max_depth = 6, eta = 0.3, subsample = 0.8, colsample_bytree = 0.8 oraz eval_metric = "logloss". Istotność cech uszeregowano według metryki gain, a następnie wybrano 30 najważniejszych genów. Las losowy zaimplementowano przy użyciu pakietu randomForest z ntree = 200. Istotność cech uszeregowano według średniego spadku współczynnika Gini i wybrano 30 najważniejszych genów. Geny wybrane przez trzy metody uczenia maszynowego poddano analizie części wspólnej, aby zidentyfikować kluczowe geny do dalszych analiz.
Budowa i ocena modelu regresji logistycznej do przewidywania ryzyka
Zbiór danych GSE54837 został losowo podzielony na zestaw treningowy (70%) i zestaw testowy (30%). Na podstawie zestawu treningowego zbudowano model regresji logistycznej przy użyciu funkcji glm z pakietu MASS, przyjmując poziomy ekspresji kluczowych genów jako cechy wejściowe. Wydajność modelu oceniono za pomocą krzywych ROC wygenerowanych przy użyciu pakietu pROC. Przedziały ufności 95% dla AUC obliczono za pomocą 2000 replik bootstrapowych. Kalibrację modelu oceniono za pomocą krzywych kalibracyjnych wygenerowanych z 1000 prób bootstrapowych (pakiet rms). Analizę DCA przeprowadzono przy użyciu pakietu dca w celu oceny wypadkowej korzyści klinicznej dla zakresu prawdopodobieństw progowych. W celu ułatwienia indywidualnej oceny ryzyka skonstruowano nomogram przy użyciu funkcji nomogram z pakietu rms.
Równanie regresji miało postać:
logit(P) = 0.5823 + 0.6010 × UPP1 - 0.6563 × PTRF + 0.3853 × B4GALT2 - 0.3972 × FAM168B + 0.1848 × PRKCDBP - 0.4787 × TOR3A. (1)
W tym przypadku P oznacza przewidywane prawdopodobieństwo wystąpienia COPD, a każdy współczynnik reprezentuje wkład odpowiadającej mu wartości ekspresji genu w logarytm szans wystąpienia COPD.
Analiza ekspresji, sieć GeneMANIA oraz molekularna sieć regulacyjna
Poziomy ekspresji genów pomiędzy grupą z POChP a grupą kontrolną w zbiorze danych GSE54837 porównano za pomocą testu sum rang Wilcoxona. W celu wizualizacji rozkładu poziomów ekspresji wygenerowano wykresy pudełkowe z wykorzystaniem pakietu ggplot2, z naniesioną medianą, rozstępem międzykwartylowym (IQR) oraz poszczególnymi punktami danych. Do konstrukcji sieci genowych i przewidywania interakcji funkcjonalnych wykorzystano narzędzie GeneMANIA. Wyszukiwanie przeprowadzono z parametrami domyślnymi: gatunek = Homo sapiens, maksymalna liczba powiązanych genów = 20. Powstałą sieć pobrano i zwizualizowano, przy czym kolory krawędzi wskazywały typy interakcji. W celu zbadania mechanizmów regulacji potranskrypcyjnej skonstruowano sieć konkurencyjnych endogennych RNA (ceRNA). miRNA celujące w sześć kluczowych genów przewidziano z wykorzystaniem dwóch niezależnych baz danych: DIANA-microT (wynik ≥ 0,8) oraz miRanda (wynik ≥ 140, energia ≤ −20 kcal/mol). Do konstrukcji par miRNA-mRNA wykorzystano część wspólną miRNA zidentyfikowanych przez obie bazy danych. Następnie, z wykorzystaniem bazy danych StarBase, przewidziano lncRNA celujące w te miRNA. Sieć regulacyjną lncRNA-miRNA-mRNA skonstruowano i zwizualizowano przy użyciu programu Cytoscape. Relacje regulacji transkrypcyjnej przewidziano za pomocą ChIP-X Enrichment Analysis Version 3 (ChEA3). Dla każdego kluczowego genu z przewidzianymi TF wybrano 10 czynników transkrypcyjnych o najwyższych wynikach wzbogacenia. Sieć regulacyjną TF-cel skonstruowano w programie Cytoscape.
Analiza wzbogacenia zestawów genów oraz ocena infiltracji komórek odpornościowych
Analizę wzbogacenia zestawów genów (GSEA) przeprowadzono przy użyciu pakietu clusterProfiler w celu zbadania funkcji biologicznych każdego kluczowego genu. Dla każdego kluczowego genu próbki podzielono na grupy o wysokiej i niskiej ekspresji w oparciu o wartość mediany. Analizę ekspresji różnicowej pomiędzy dwiema grupami wykonano za pomocą pakietu limma, a otrzymaną listę genów uszeregowano według wartości signed log₂ fold-change. Analizę GSEA przeprowadzono przy użyciu funkcji gseGO dla terminów procesów biologicznych GO oraz funkcji gseKEGG dla szlaków KEGG, stosując następujące parametry: minGSSize = 10, maxGSSize = 500, pvalueCutoff = 0.05 oraz nPerm = 1 000. Względną obfitość 28 typów komórek odpornościowych oszacowano za pomocą analizy wzbogacenia zestawów genów dla pojedynczych próbek (ssGSEA) zaimplementowanej w pakiecie GSVA. Kuratowana macierz sygnatur zestawów genów zawierająca geny markerowe dla 28 typów komórek odpornościowych została pozyskana z literatury przedmiotowej19. Dla każdej próbki zastosowano funkcję gsva z parametrami method = "ssgsea", ssgsea.norm = TRUE oraz kcdf = "Gaussian". Współczynniki korelacji Spearmana pomiędzy wynikami wzbogacenia ssGSEA a poziomami ekspresji sześciu kluczowych genów obliczono za pomocą funkcji cor.test. Wartości p skorygowano pod kątem wielokrotnych testów metodą Benjamini-Hochberga. Macierz korelacji zwizualizowano w formie mapy ciepła (heatmap) przy użyciu pakietu pheatmap.
Przewidywanie leków, dokowanie molekularne i analiza powiązań z chorobami
Potencjalne związki terapeutyczne ukierunkowane na kluczowe geny zidentyfikowano przy użyciu bazy danych DrugBank. W programie Cytoscape skonstruowano sieć interakcji „lek ukierunkowany na kluczowy gen”, aby zwizualizować przewidywane oddziaływania lek-gen. Przeprowadzono dokowanie molekularne przy użyciu platformy CB-Dock2 w celu oceny powinowactwa wiązania. Trójwymiarową strukturę białka ludzkiego UPP1 pobrano z Protein Data Bank (PDB ID: 7B8T). Struktury molekularne leków (w formacie SMILES) pozyskano z bazy PubChem. Dokowanie wykonano przy użyciu silnika AutoDock Vina, a wyniki uszeregowano według wolnej energii wiązania (ΔG, w kcal/mol). Kompleksy dokowania zwizualizowano za pomocą programu PyMOL. Powiązania między kluczowymi genami a chorobami ludzi związanymi z ekspozycją na czynniki środowiskowe zbadano przy użyciu Comparative Toxicogenomics Database (CTD). Każdy gen analizowano oddzielnie, a następnie wyekstrahowano dziesięć najsilniej powiązanych chorób i przedstawiono je za pomocą wykresów radarowych.
Protokół RT-qPCR
Próbki krwi żylnej obwodowej pobrano od ośmiu pacjentów z POChP oraz ośmiu zdrowych osób z grupy kontrolnej w szpitalu Shenzhen Luohu Hospital of Traditional Chinese Medicine. POChP zdiagnozowano zgodnie z kryteriami Global Initiative for Chronic Obstructive Lung Disease (GOLD), definiującymi chorobę jako stosunek FEV1/FVC < 0,70 po podaniu leku rozszerzającego oskrzela. Grupa kontrolna składała się z dopasowanych pod względem wieku i płci zdrowych ochotników, bez historii chorób dróg oddechowych i z prawidłowymi wynikami prób czynnościowych płuc (przewidywane FEV1% ≥ 80% oraz FEV1/FVC ≥ 0,70). Informacje wyjściowe dotyczące pacjentów przedstawiono w Tabeli 2. Całkowite RNA wyekstrahowano z próbek krwi pacjentów z POChP przy użyciu zestawu do ekstrakcji RNA z krwi. W celu syntezy cDNA, 500 ng całkowitego RNA poddano odwrotnej transkrypcji z wykorzystaniem zestawu do syntezy cDNA z usuwaniem genomowego DNA, zgodnie z dostarczonym protokołem. Otrzymane cDNA rozcieńczono do stężenia 150 ng/μL.
Analizę RT-qPCR przeprowadzono z wykorzystaniem qPCR master mix opartego na SYBR Green w systemie PCR w czasie rzeczywistym. Każda reakcja o objętości 10 μL zawierała 5 μL 2x SYBR Green master mix, po 0,5 μL starterów bezpośrednich i odwrotnych (10 μM), 1 μL rozcieńczonego cDNA (15 ng/μL) oraz 3 μL wody wolnej od nukleaz. Warunki cykliczności obejmowały wstępną denaturację w 95 °C przez 5 min, a następnie 40 cykli: 95 °C przez 10 s i 60 °C przez 30 s, z końcową analizą krzywej topnienia od 60 °C do 95 °C w celu weryfikacji swoistości amplifikacji. Wszystkie reakcje przeprowadzono w trzech powtórzeniach technicznych. Jako wewnętrzny gen referencyjny wykorzystano β-actin. Efektywność starterów dla każdego genu docelowego została zwalidowana za pomocą serii rozcieńczeń krzywej standardowej i mieściła się w zakresie od 90% do 110%. Poziomy ekspresji genów znormalizowano do β-actin, a względną ekspresję obliczono metodą 2-ΔΔCt. Porównania statystyczne między grupą z COPD a grupą kontrolną przeprowadzono za pomocą testu U Manna-Whitneya.
Analiza statystyczna
Wizualizacje sieci zostały utworzone przy użyciu programu Cytoscape, a analizy statystyczne przeprowadzono w oprogramowaniu R. O ile nie wskazano inaczej, w przypadku danych o rozkładzie nienormalnym zastosowano test U Manna-Whitneya, a w przypadku danych o rozkładzie normalnym test t Studenta w celu porównania dwóch grup. Wartość p < 0,05 uznano za statystycznie istotną.