Artykuł metodologiczny
* These authors contributed equally
Ten protokół przedstawia standaryzowany przepływ pracy bioinformatycznej do analizy zmian transkryptomowych w ostrym białaczu szpiaku (AML). Celem jest porównanie nowo rozpoznanych i nawrotowych próbek szpiku kostnego oraz uporządkowanie pod względem priorytetu sygnatur molekularnych związanych z opornością na chemioterapię i progresją choroby do dalszego badania.
Ostra białaczka szpikowa (AML) jest wysoce heterogeniczną chorobą hematologiczną, w której nawrót i nabyte oporność na chemioterapię pozostają głównymi przyczynami niepowodzenia leczenia. W tym artykule przedstawiono protokół bioinformatyczny do analizy transkryptomowej aspirowanych próbek szpiku kostnego. Głównym celem protokołu jest zapewnienie standaryzowanego przepływu pracy w celu zidentyfikowania sygnatur molekularnych związanych z progresją choroby i opornością na terapię w nawróconej AML. Przepływ opisuje procedury obliczeniowe do porównywania niesparowanych próbek szpiku kostnego, demonstrując je na danych sekwencjonowania z pięciu nowo rozpoznanych przypadków i czterech przypadków nawrotu. Ta metoda opisuje kluczowe kroki przetwarzania danych sekwencjonowania RNA, przeprowadzania analizy różnicowej ekspresji genów i przeprowadzania dalszych ewaluacji funkcjonalnych. Zastosowanie tego przepływu pozwoliło na zidentyfikowanie 2 025 różnicowej ekspresji genów (DEGs), w tym FOXC1, HOXA11, HOXA11-AS i AXL, jako kandydatów do transkryptów związanych z nawrotami w tym reprezentatywnym zestawie danych. Analizy funkcjonalne i sieciowe dodatkowo uporządkowały zbiory genów i centra interakcji związane z sygnałowaniem małych GTPaz, sygnałowaniem zapalnym, interakcjami z macierzą pozakomórkową i procesami biosyntezy RNA. Ogólnie rzecz biorąc, ta metodologia zapewnia powtarzalny komputerowy przepływ do mapowania sygnatur transkryptomowych związanych z nawróconej AML i generowania hipotez wymagających dalszej eksperymentalnej walidacji.
Ostra białaczka szpikowa (AML) to grupa klonalnych złośliwych nowotworów wywodzących się z komórek macierzystych/przyprostokształtnych układu krwiotwórczego, charakteryzująca się nieprawidłową proliferacją niedojrzałych komórek szpikowych w rdzeniu kostnym i tłumieniem różnicowania układu krwiotwórczego1,2,3. Mimo że obecnie standardowa chemioterapia indukcyjna (jak cytarabina połączona z antracyklinami) może wywołać pełną remisję u większości pacjentów, wskaźnik nawrotów pozostaje wysoki, wynosząc 50–70%, a rokowanie dla pacjentów z nawrotami jest znacznie gorsze4,5,6. Nabyta oporność na chemioterapię jest związana z niepowodzeniem leczenia AML, co wymaga głębokiej analizy sygnałów molekularnych związanych z tym procesem, aby poprawić strategie terapeutyczne i wskaźniki przeżycia pacjentów7,8.
W szerszym korpusie literatury, istniejące badania wskazują, że oporność na chemioterapię w AML nie ogranicza się do nadregulacji pomp wydalających leki lub nieprawidłowego metabolizmu leków, ale jest również związana z utrzymaniem przeżycia komórek macierzystych białaczki (LSCs), tworzeniem fenotypów przypominających przemianę nabłonkowo-mezenchymalną (EMT) w środowisku hematologicznym oraz interakcjami z mikrośrodowiskiem rdzenia kostnego9,10,11. Na przykład, populacje LSC wykazują wysoką zdolność do samoregeneracji i stan spoczynku, co wiąże się z wrodzoną odpornością na środki chemioterapeutyczne działające w określonych fazach cyklu komórkowego12. Ponadto, nadregulacja receptorów kinaz tyrozynowych, takich jak AXL, została związana z opornością w przypadku AML z mutacją FLT3-ITD, występującą wraz z aktywacją szlaków PI3K/AKT i MAPK oraz zwiększoną zdolnością anty-apoptotyczną13,14. Przeprogramowanie metaboliczne i remodelacja epigenetyczna zostały również rozpoznane jako ważne osie regulacyjne w formowaniu oporności. Dowody sugerują, że komórki AML podczas nawrotu mogą przystosowywać się do stresu oksydacyjnego i uszkodzenia DNA wywołanego chemioterapią poprzez zwiększoną aktywność fosforylacji oksydacyjnej (OXPHOS), modulowane proporcje NAD⁺/NADH oraz zmienione stany modyfikacji histonów15,16,17. Czynniki zapalne w mikrośrodowisku rdzenia kostnego, takie jak IL-6 i CXCL8, są również związane z przeżyciem LSC i opornością na chemioterapię, często występujące w koordynacji z aktywacją szlaków sygnalizacyjnych STAT3/NF-κB11,18.
Mimo tych rozpoznanych mechanizmów, zmiany transkryptomowe związane z nawrotami AML i opornością na chemioterapię pozostają niedostatecznie scharakteryzowane, zwłaszcza jeśli oceniana jest łącznie metabolizm, epigenetyka i sygnały związane z mikrośrodowiskiem rdzenia kostnego w próbkach klinicznych. Ten przepływ pracy spełnia tę potrzebę, priorytetyzując kandydatów do genów różniczkowo ekspresowanych (DEGs), szlaków i sieci regulacyjnych związanych z przejściem od początkowej diagnozy do klinicznego nawrotu. Podejście to integruje analizę różnicowej ekspresji, analizę wzbogacania zestawów genów (GSEA) oraz budowę sieci interakcji białko-białko (PPI) w celu odwzorowania przeprogramowania transkrypcyjnego na skalę systemową i wskazania kandydatów do dalszej walidacji mechanotycznej.
Ogólnym celem tej metody jest przedstawienie standaryzowanego, powtarzalnego pipeline'u bioinformatycznego do porównywania transkryptomów masowych nowo rozpoznanych i nawrotowych próbek szpiku z AML. Racjonalnym zastosowania tej techniki in silico jest jej zdolność do rejestracji bezstronne, genomowe wydarzenia transkrypcyjne, wykraczając poza ograniczenia analiz pojedynczych szlaków, aby priorytetyzować złożone, wielowymiarowe sieci regulacyjne. Ta technika oferuje znaczne zalety w porównaniu z alternatywnymi metodami, takimi jak mikroarrays czy ukierunkowane panele qPCR wielokrotnego zmapu, poprzez zapewnienie wyższego zakresu dynamicznego, zdolności do wykrywania nowych transkryptów oraz precyzyjnego ilorazowanego ekspresji genów bez ograniczeń związanych z wstępnie zaprojektowanymi sondami19,20. Aby ustalić, czy ta metoda jest odpowiednia dla ich zastosowania, czytelnicy powinni zauważyć, że ten pipeline jest specjalnie zaprojektowany dla badaczy przetwarzających dane z sekwencjonowania RNA w ilości masowej z sparowanych lub niesparowanych kohortach klinicznych, takich jak aspiraty tkankowe. Jest odpowiedni do identyfikacji szerokich sygnatur związanych z opornością i kandydatów do sieci regulacyjnych, podczas gdy naukowcy wymagający rozdzielczości specyficznej dla typu komórki lub przestrzennej muszą skorzystać z uzupełniających przepływów sekwencjonowania jednokomórkowego lub przestrzennego. O
Wszystkie metody wykorzystujące pobieranie próbek tkanki ludzkiej zostały wykonane zgodnie z wytycznymi instytucjonalnymi i Deklaracją Helsińską (zmieniona w 2013 roku). Kliniczne próbki szpiku kostnego zostały uzyskane za zgodą instytucjonalnego komitetu etycznego (Numery zgody TY-ZKY2024-116-01 i TY-ZKY2024-116-02).
1. Pobieranie próbek klinicznych i klasyfikacja pacjentów
2. Kontrola jakości RNA i przygotowanie biblioteki
Kliniczne badanie pacjentów i walidacja sekwencjonowania
Pomyślne przeprowadzenie protokołu ekstrakcji RNA w górnym strumieniu i przygotowania bibliotek (Rysunek 1) zostało potwierdzone przez wydajność sekwencjonowania i metryki jakości. W tym reprezentatywnym zestawie danych próbki szpiku kostnego od pięciu nowo rozpoznanych pacjentów z AML i czterech pacjentów z nawrotowym AML dały średnio około 6,0 GB surowych danych na próbkę. Ocena kontroli jakości (Tabela 2) potwierdziła, że jakość zasad i głębokość odczytu spełniały progi wymagane do dalszej analizy bioinformatycznej9. Niska jakość RNA (na przykład RIN < 6,0), niska stopa mapowania lub wysokie zniekształcenie transkryptu byłyby oznaką suboptymalnej jakości wejścia i kompromitująby niezawodność dalszej analizy różnicowej ekspresji.
Globalna wariancja transkryptomu i PCA
Aby ocenić globalną wariancję transkryptomu i sprawdzić klasteryzację kliniczna, PCA została przeprowadzona na znormalizowanych danych ekspresji. W tym reprezentatywnym zestawie danych, grupy nowo rozpoznanych i nawrotowych wykazywały rozdzielenie w przestrzeni dwuwymiarowej (Rysunek 2A)20, przy czym PC1 i PC2 odpowiadały za 23,82% i 18,75% całkowitej wariancji, odpowiednio. Diagramy Venna w Rysunku 2B,C zapewniają dodatkowe opisowe podsumowanie genów wykrytych w próbkach w grupach nowo rozpoznanych i nawrotowych, wspierające sprawdzanie powtarzalności na poziomie próbek przed dalszą analizą różnicowej ekspresji. Ponieważ kohorta była niewielka i nieparzysta, separacja PCA została przedstawiona jako ilustracyjne wyjście przepływu pracy, a nie jako ostateczny dowód specyficznej dla stanu choroby biologii.
Analiza genów o zróżnicowanej ekspresji (DEG)
Zastosowanie ustalonych progów protokołu (|log2FC| ≥ 1 i wartość P dostosowana ≤ 0,05) na wyjściu DESeq2 pozwoliło na zidentyfikowanie 2025 genów DEG, w tym 772 genów nadekspresowanych i 1253 genów podekspresowanych w grupie nawrotowej (Rysunek 3A). Kandydaci do transkryptów o dużej zmienności obejmowali FOXC1 (log2FC = 7,55, P = 4,92 x 10-5), HOXA11 (log2FC = 7,76), HOXA11-AS (log2FC = 7,23) i AXL (log2FC = 3,50), wraz z podekspresją RHOB, PTX3 i CXCL8. Istniejąca literatura łączy kilka z tych genów z komórkami macierzystymi AML, sygnalizacją lub odpowiedzią na terapię13,29; jednak obecny przepływ pracy identyfikuje je jedynie jako kandydaci do transkryptów związanych z nawrotami. Jakakolwiek jednoznaczny mechanizm roli w oporności klinicznej wymaga późniejszej niezależnej weryfikacji funkcjonalnej.
Funkcjonalne i ścieżkowe wzbogacenie (GO, KEGG i GSEA)
Protokół analizy funkcjonalnej zmapowywał geny DEG na szersze systemy biologiczne. Analiza GO wykryła wzbogacenie terminów związanych z sygnalizacją z pośrednictwem małych GTPaz, transportem jonów metali i składaniem chromatyny (Rysunek 4A-C). Mapowanie ścieżek KEGG wykryło skojarzenia z interakcjami receptorów ECM i interakcjami cytokin-cytokinowych receptorów (Rysunek 4D). GSEA wykazała wzbogacenie procesów biosyntezy RNA w grupie nawrotowej i wzbogacenie ścieżek metabolizmu energetycznego w grupie nowo rozpoznanych (Rysunek 5A). Te wyniki wzbogacenia dostarczają opisową mapę zmienionych zestawów genów i powinny być interpretowane jako skojarzenia generujące hipotezy, a nie udowodnione czynniki odpowiedzialne za nawroty.
Budowa sieci interakcji białko-białko (PPI)
Początkowa sieć STRING zawierała 56 węzłów i 193 interakcje. Po usunięciu odłączonych lub osieroconych węzłów, wyświetlona pod-sieć Cytoscape zawierała 42 węzły i 136 interakcji (Rysunek 5B). Analiza modularna sieci priorytetowała TP53, CCL2, CXCL8 i IL6 jako centralne matematyczne węzły z największą liczbą interakcji. Ponieważ sieć PPI opiera się na przewidywanych wynikach interakcji z bazy danych (np. wynik ATF3: 0,982), identyfikacja węzłów powinna być interpretowana jako priorytetyzacja celów dla przyszłych badań empirycznych, a nie jako bezpośredni dowód uniknięcia apoptozy przez p53 lub innych mechanizmów oporu.
Dane surowe z sekwencjonowania RNA wygenerowane w tym protokole zostały zdeponowane w repozytorium Figshare i są publicznie dostępne za pośrednictwem następującego DOI: https://doi.org/10.6084/m9.figshare.30655814. Przetworzone dane i powiązane pliki analityczne są zawarte w artykule i jego materiałach uzupełniających. Reprezentatywne parametry wiersza poleceń i ustawienia analityczne użyte do odtworzenia przepływu pracy obliczeniowego są dostarczane jako Plik uzupełniający 1.
Kluczowe kroki w protokole
Pomyślne wykonanie tego przepływu pracy bioinformatycznej opiera się na kilku kluczowych krokach. Po pierwsze, natychmiastowe zamrażanie i odpowiednie lizysowanie aspyratu szpiku kostnego (Krok 1.6) są zasadnicze, ponieważ mikrośrodowisko szpiku kostnego jest bogate w rybonukleazy, które mogą szybko niszczyć integralność transkryptomu30. Podczas fazy obliczeniowej, poprawna konfiguracja wzorca projektu eksperymentu w pakiecie DESeq2 (Krok 6.3) jest krytyczna dla dokładnego różnicowania ekspresji, szczególnie gdy porównywany jest stan kliniczny (nowo rozpoznany w porównaniu z nawracanym), podczas gdy brane są pod uwagę potencjalne zmienne mylące. Wreszcie, zastosowanie rygorystycznych progów współczynnika fałszywych odkryć (FDR) podczas Analizy Obogacenia Zestawów Genów (GSEA) (Krok 9.6) jest kluczowy statystycznym kontrolnym, aby zapobiec nadmiernej interpretacji fałszywych dodatnich funkcjonalnych sieci.
Modyfikacje i rozwiązywanie problemów
Częstym wyzwaniem w tej metodzie jest obecność efektów partii, które często występują, gdy próbki kliniczne są zbierane i sekwencjonowane w dłuższych przedziałach czasowych. Zmienne partii powinny być zdefiniowane przed analizą, w tym data pobrania próbki, data ekstrakcji RNA, partia przygotowania biblioteki, pasmo sekwencjonowania i uruchomienie sekwencjonowania. Jeśli analiza PCA lub korelacji próbek ujawnia grupowanie na podstawie daty sekwencjonowania lub innej zmiennej technicznej, a nie na podstawie fenotypu klinicznego, użytkownicy powinni zmodyfikować protokół, włączając zmienną partii do wzorca różnicowej ekspresji, gdy jest to statystycznie wykonalne lub zastosować algorytmy korekcji partii, takie jak ComBat lub SVA, przed wizualizacja31. Jeśli stosuje się ten protokół do pełnej krwi zamiast aspyratów szpiku kostnego, istotną modyfikacją jest włączenie kroków wyczerpania mRNA globiny podczas przygotowania biblioteki, aby zapobiec monopolizacji głębokości odczytu sekwencjonowania przez wysoce obfite transkrypty globiny. Wersje oprogramowania i główne parametry reprezentatywnego przepływu pracy zostały uzupełnione następująco: fastp v0.23.2, HISAT2 v2.0.5, StringTie v1.3.3b, featureCounts v1.5.0-p3, R v3.5.0, DESeq2 v1.20.0, clusterProfiler v3.8.1, org.Hs.eg.db v3.6.0, sekwencjonowanie paired-end 150 bp, GSEA v4.2.3 z 1 000 permutacjami zestawów genów, MSigDB v7.5.1, STRING v11.5 z interakcjami o wysokiej pewności, i Cytoscape v3.9.1. Reprezentatywne parametry wiersza poleceń i ustawienia analizy są dostępne w Pliku uzupełniającym 1.
Ograniczenia metody
Mimo że jest kompleksowy, ten protokół ma wrodzone ograniczenia metodologiczne. Po pierwsze, wykorzystuje sekwencjonowanie RNA masowego, którego uchwytuje średni profil transkryptomu całego aspyratu szpiku kostnego i brakuje rozdzielczości przestrzennej pojedynczych komórek. Dlatego przepływ pracy nie może określić, czy wzmożony podpis związany z nawracaniem pochodzi od komórek macierzystych białaczki, komórek stromalnych, komórek odpornościowych, czy też zmian w składzie typu komórki32. Po drugie, reprezentatywny zestaw danych jest mały (n = 9) i niesparowany, co ogranicza moc statystyczną i uniemożliwia jednoznaczne wyciąganie wniosków przyczynowych. Po trzecie, przepływ pracy jest czysto in silico. Generuje kandydatów do regulatorowych węzłów i ścieżek sygnałowych, ale nie może niezależnie zweryfikować ich funkcjonalnej konieczności w chemorezystancji bez ortogonalnej walidacji eksperymentalnej in vitro lub in vivo.
Ostatnie badania pojedynczych komórek i genomów pojedynczych komórek poszerzyły ramy odniesienia AML, rozdzielając heterogeniczność stanu komórki, architekturę klonalną i ewolucję związaną z terapią z wyższą rozdzielczością33,34,35,36. Te podejścia uzupełniają opisany tu przepływ pracy sekwencjonowania RNA masowego: sekwencjonowanie masowe zapewnia praktyczną i opłacalną strategię przesiewową dla transkryptomu na poziomie kohorty, podczas gdy metody pojedynczych komórek i multi-omiczne mogą być stosowane w badaniach następczych w celu przypisania kandydatów do sygnałów do określonych populacji złośliwych lub mikrośrodowiskowych komórek.
Znaczenie w stosunku do istniejących metod
Mimo tych ograniczeń, ten przepływ pracy transkryptomu oferuje zalety nad alternatywnymi technikami diagnostycznymi i analitycznymi. Tradycyjne oceny kliniczne nawrotu AML często polegają na docelowych panelach qPCR multiplex lub standardowej cytometrii przepływowej. Chociaż są przydatne do szybkiego rozpoznawania, te ukierunkowane metody są ograniczone przez wstępnie zdefiniowane sondy i mogą oceniać tylko znane markery oporu19. Wykorzystując bezstronne sekwencjonowanie transkryptomu na skalę genomu połączone z analizą sieci, ten protokół może nominować nowe transkrypty i asocjacje
Autorzy deklarują brak konfliktu interesów.
Badanie to było sfinansowane przez Biuro Nauki i Technologii Miasta Ganzhou (2022—ZD1368).
| Nazwa | Firma | Numer katalogowy | Komentarze |
|---|---|---|---|
| Agilent 2100 Bioanalyzer | Agilent Technologies, Santa Clara, CA, USA | RRID:SCR_019389 | G2939BA |
| AMPure XP system | Beckman Coulter, Brea, CA, USA | RRID:SCR_008452 | A63881 |
| cBot Cluster Generation System | Illumina, San Diego, CA, USA | - | SY-301-2002 lub identyfikator systemu specyficznego dla instytucji |
| clusterProfiler (oprogramowanie) | Bioconductor | RRID:SCR_016884 | v3.8.1 |
| Cytoscape (oprogramowanie) | Cytoscape Consortium | RRID:SCR_003032 | v3.9.1 |
| DESeq2 (oprogramowanie) | Bioconductor | RRID:SCR_015687 | v1.20.0 |
| DNA Polymerase I | New England Biolabs (NEB, Ipswich, MA, USA) | - | M0209L |
| dNTP Solution Mix | New England Biolabs (NEB, Ipswich, MA, USA) | - | N0447L |
| edgeR (oprogramowanie) | Bioconductor | RRID:SCR_012802 | v3.22.5 |
| fastp (oprogramowanie) | OpenGene | RRID:SCR_016962 | v0.23.2 |
| featureCounts / Subread (oprogramowanie) | The Walter and Eliza Hall Institute | RRID:SCR_012919 | featureCounts v1.5.0-p3 |
| GRCh38 reference genome | Genome Reference Consortium / Ensembl | - | GRCh38; Ensembl release 109 |
| GSEA software | Broad Institute | RRID:SCR_003199 | v4.2.3 |
| HISAT2 (oprogramowanie) | Johns Hopkins University | RRID:SCR_015530 | v2.0.5 |
| M-MuLV Reverse Transcriptase (RNase H-) | New England Biolabs (NEB, Ipswich, MA, USA) | - | M0253L |
| MSigDB gene sets | Broad Institute | RRID:SCR_016863 | v7.5.1 |
| NEBNext Ultra II Directional RNA Library Prep Kit for Illumina | New England Biolabs (NEB, Ipswich, MA, USA) | - | E7760L/E7765L lub zestaw specyficzny dla laboratorium |
| NEBNext Ultra II RNA Library Prep Kit for Illumina | New England Biolabs (NEB, Ipswich, MA, USA) | - | E7770L |
| NovaSeq sequencing platform | Illumina, San Diego, CA, USA | RRID:SCR_016387 | NovaSeq system; identyfikator instrumentu dostawcy usług |
| org.Hs.eg.db (pakiet adnotacji) | Bioconductor | - | v3.6.0 |
| Phusion High-Fidelity DNA Polymerase | Thermo Fisher Scientific, Waltham, MA, USA | RRID:AB_2756816 | F530L |
| Qubit 2.0 Fluorometer | Thermo Fisher Scientific, Waltham, MA, USA | RRID:SCR_018095 | Q32866 |
| Qubit dsDNA HS Assay Kit | Thermo Fisher Scientific, Waltham, MA, USA | - | Q32851 |
| R software | R Foundation for Statistical Computing | RRID:SCR_001905 | v3.5.0 |
| Random Hexamer Primer | Thermo Fisher Scientific, Waltham, MA, USA | - | SO142 |
| RNA 6000 Nano Kit | Agilent Technologies, Santa Clara, CA, USA | - | 5067-1511 |
| RNase H | New England Biolabs (NEB, Ipswich, MA, USA) | - | M0297L |
| RNA-seq Library Prep Kit / Sequencing Service | Novogene, Beijing, China | - | Nr projektu X101SC25054246-Z01-J003 |
| STRING database | STRING Consortium | RRID:SCR_005223 | v11.5 |
| StringTie (oprogramowanie) | Johns Hopkins University / Center for Computational Biology | RRID:SCR_016323 | v1.3.3b |
| TruSeq PE Cluster Kit v3-cBot-HS | Illumina, San Diego, CA, USA | - | PE-401-3001 |
Poproś o pozwolenie na ponowne wykorzystanie tekstu lub ilustracji tego artykułu JoVE
Poproś o pozwolenie