Transcriptomische profilering en bio-informatische analyse van beenmergmonsters om chemotherapieresistentiesignaturen bij acute myeloïde leukemie te identificeren

Samenvatting

Dit protocol biedt een gestandaardiseerde bio-informatische workflow om transcriptomische veranderingen bij acute myeloïde leukemie (AML) te analyseren. Het doel is om nieuw gediagnosticeerde en teruggevallen beenmergmonsters te vergelijken en moleculaire signaturen die geassocieerd zijn met chemotherapieresistentie en ziekteprogressie te prioriteren voor vervolgonderzoek.

Samenvatting

Acute myeloïde leukemie (AML) is een sterk heterogene hematologische maligniteit waarbij terugval en verworven chemoresistentie de belangrijkste oorzaken van het falen van de behandeling blijven. Dit artikel presenteert een bio-informatisch protocol voor transcriptomische analyse van beenmergaspiraten. Het primaire doel van het protocol is het bieden van een gestandaardiseerde workflow om moleculaire signaturen te identificeren die verband houden met ziekteprogressie en therapieresistentie bij terugvallende AML. De pijplijn beschrijft de computationele procedures voor het vergelijken van ongepaarde beenmergmonsters, aangetoond met sequencinggegevens van vijf nieuw gediagnosticeerde gevallen en vier terugvallende gevallen. Deze methode beschrijft de essentiële stappen voor het verwerken van RNA-sequencinggegevens, het uitvoeren van differentiële genexpressie-analyse en het uitvoeren van latere functionele evaluaties. Door deze workflow toe te passen werden 2.025 differentieel expressieve genen (DEGs), waaronder FOXC1, HOXA11, HOXA11-AS en AXL, geïdentificeerd als kandidaat-transcripten geassocieerd met terugval in deze representatieve dataset. Functionele en netwerkanalyses gaven verder prioriteit aan gensets en interactiehubs die verband houden met kleine GTPase-signalering, inflammatoire signalering, extracellulaire matrixinteracties en RNA-biosynthetische processen. Al met al biedt deze methodologie een reproduceerbare computationele pijplijn voor het in kaart brengen van transcriptomische signaturen die samenhangen met terugvallende AML en voor het genereren van hypothesen die latere experimentele validatie vereisen.

Inleiding

Acute myeloïde leukemie (AML) is een groep klonale maligne neoplasma's afkomstig van hematopoëtische stam-/voorlopercellen, gekenmerkt door abnormale proliferatie van onrijpe myeloïde cellen in het beenmerg en onderdrukking van hematopoëtische differentiatie ^1,2,3. Hoewel de huidige standaard inductiechemotherapie (zoals cytarabine gecombineerd met anthracyclines) bij de meeste patiënten volledige remissie kan induceren, blijft het terugvalpercentage zo hoog als 50%–70%, en is de prognose voor teruggevallen patiënten significant verslechterd ^4,5,6. Verworven chemotherapieresistentie wordt geassocieerd met het falen van de AML-behandeling, waardoor een diepgaande analyse van de moleculaire signaturen die aan dit proces gekoppeld zijn noodzakelijk is om therapeutische strategieën en patiëntoverlevingskansen te verbeteren ^7,8.

In de bredere literatuur geven bestaande studies aan dat chemotherapieresistentie bij AML niet beperkt is tot de upregulatie van effluxpompen of abnormale medicijnmetabolisme, maar ook verband houdt met het behoud van leukemie-stamcellen (LSC's), de vorming van epitheel-mesenchymale overgangsfenotypes (EMT) binnen de hematologische niche, en interacties met het beenmergmicro-milieu^{9, 10,11}. Zo vertonen LSC-populaties een hoge zelfvernieuwingscapaciteit en ruststand, wat geassocieerd is met inherente resistentie tegen celcyclus-specifieke chemotherapeutische middelen¹². Daarnaast is de upregulatie van receptor-tyrosinekinasen, zoals AXL, geassocieerd met resistentie in FLT3-ITD+ AML, die plaatsvindt naast de activatie van PI3K/AKT en MAPK-routes en verbeterde anti-apoptotische mogelijkheden^13,14. Metabole herprogrammering en epigenetische remodellering zijn ook erkend als belangrijke regulerende assen bij de vorming van resistentie. Bewijs suggereert dat AML-cellen tijdens terugval zich kunnen aanpassen aan door chemotherapie veroorzaakte oxidatieve stress en DNA-schade door verhoogde oxidatieve fosforylering (OXPHOS)-activiteit, gemoduleerde NAD⁺/NADH-verhoudingen en veranderde histonmodificatietoestanden 15,16,17. Ontstekingsfactoren in het beenmergmicro-omgeving, zoals IL-6 en CXCL8, zijn ook geassocieerd met LSC-overleving en chemotherapieresistentie, vaak in samenwerking met de activatie van STAT3/NF-κB signaalroutes^11,18.

Ondanks deze erkende mechanismen blijven de transcriptomische veranderingen die geassocieerd worden met AML-terugval en chemoresistentie onvolledig gekarakteriseerd, vooral wanneer metabole, epigenetische en aan het beenmerg gerelateerde micro-omgevingssignalen samen in klinische monsters worden geëvalueerd. Deze workflow pakt deze behoefte aan door kandidaat-DEGs, paden en regulerende netwerken te prioriteren die samenhangen met de overgang van de eerste diagnose naar klinische terugval. De benadering integreert differentiële expressieanalyse, gensetverrijkingsanalyse (GSEA) en eiwit-eiwitinteractie (PPI) netwerkconstructie om systeembrede transcriptieherprogrammering in kaart te brengen en kandidaten te nomineren voor latere mechanistische validatie.

Het algemene doel van deze methode is het presenteren van een gestandaardiseerde, reproduceerbare bio-informatica-pijplijn voor het vergelijken van de bulk transcriptomen van nieuw gediagnosticeerde versus recidiverende AML-beenmergmonsters. De reden voor het gebruik van deze in silico-techniek is het vermogen om onbevooroordeelde, genoombrede transcriptiegebeurtenissen vast te leggen, waarmee de beperkingen van enkelvoudige padanalyses worden overschreden om complexe, multidimensionale regulerende netwerken te prioriteren. Deze techniek biedt aanzienlijke voordelen ten opzichte van alternatieve methoden, zoals microarrays of gerichte multiplex qPCR-panelen, door een hoger dynamisch bereik te bieden, de mogelijkheid om nieuwe transcripten te detecteren en precieze kwantificering van genexpressie zonder de beperkingen van vooraf ontworpen probes^19,20. Om te bepalen of deze methode geschikt is voor hun toepassing, moeten lezers opmerken dat deze pijplijn specifiek is ontworpen voor onderzoekers die bulk RNA-sequencinggegevens verwerken van gepaarde of ongepaarde klinische cohorten, zoals weefselaspiraten. Het is geschikt voor het identificeren van brede weerstandsgerelateerde signaturen en kandidaat-regulatorische netwerken, terwijl onderzoekers die celtype-specifieke of ruimtelijke resolutie nodig hebben, complementaire single-cell of ruimtelijke sequencing-workflows moeten gebruiken. Uiteindelijk maakt dit computationele protocol het mogelijk om kandidaatgenen, routes en regulerende netwerken te prioriteren voor latere experimentele onderzoeken.

Protocol

Alle methoden waarbij menselijk weefselmonster wordt gebruikt, zijn uitgevoerd in overeenstemming met institutionele richtlijnen en de Verklaring van Helsinki (herzien in 2013). De klinische beenmergmonsters werden verkregen met goedkeuring van de ethische commissie van de instelling (goedkeuringsnummers TY-ZKY2024-116-01 en TY-ZKY2024-116-02).

1. Klinische monsterverzameling en patiëntclassificatie

1. Selecteer beenmergaspiratiemonsters van patiënten die formeel zijn gediagnosticeerd met acute myeloïde leukemie (AML) op basis van de classificatiecriteria van de Wereldgezondheidsorganisatie (WHO).
2. Pas specifieke inclusie- en exclusiecriteria toe tijdens de patiëntselectie om cohorte-homogeniteit en reproduceerbaarheid te waarborgen. Neem volwassen patiënten met primaire AML op en sluit patiënten met secundaire AML, acute promyelocytaire leukemie of een eerdere voorgeschiedenis van andere maligniteiten uit (Tabel 1).
3. Wijs het verzamelde monster toe aan de nieuw gediagnosticeerde groep als de patiënt bij de eerste klinische diagnose onbehandelde AML vertoont.
4. Wijs het verzamelde monster toe aan de terugvallende groep als de patiënt een terugkeer van leukemische blasts in het perifere bloed of meer dan 5% blasts in het beenmerg vertoont na een gedocumenteerde volledige remissie.
5. Verzamel de gedeïdentificeerde resterende beenmergmonsters direct na de routinematige klinische beenmergaspiratieprocedure.
   OPMERKING: Bij de demonstratie van dit specifieke protocol werden van september 2024 tot september 2025 negen opeenvolgende exemplaren (vijf nieuw gediagnosticeerd en vier terugvallen) verzameld. Omdat hiermee alleen gedeïdenificeerde resterende klinische monsters werden gebruikt, schafte de ethische commissie de verplichting van schriftelijke geïnformeerde toestemming op.
6. Verwerk onmiddellijk de verzamelde beenmergaspiratie voor RNA-behoud met behulp van de standaard fenol-guanidinium-lysemethode²¹.
   1. Breng het vers verzamelde beenmergaspirat over in een verzamelbuis met een antistollingsmiddel. Schud de buis om de aspiratie en het antistollingsmiddel grondig te mengen.
   2. Extraheren van een gemeten volume van het anticoaguleerde monster en voeg dit toe aan een commercieel fenol-guanidinium-lysereagens. Houd een volumeverhouding van 3 delen lysisreagens op 1 deel monster.
      VOORZICHTIG: Het fenol-guanidinium lyse-reagens bevat giftige en corrosieve chemicaliën die ernstige brandwonden en weefselschade kunnen veroorzaken. Voer alle reagentia-omgang uit in een chemische dampkap terwijl je de juiste persoonlijke beschermingsmiddelen draagt.
   3. Schud de buis krachtig om het monster en het lysisreagens volledig te homogeniseren. Zorg dat het mengsel volledig is gemengd en controleer dat er geen zichtbare stolsels in de oplossing achterblijven.
   4. Vries het gehomogeniseerde monster onmiddellijk in een snap-frees door de buis onder te dompelen in vloeibare stikstof.
      VOORZICHTIG: vloeibare stikstof is extreem koud en kan ernstige bevriezing veroorzaken bij contact. Draag cryogene handschoenen en een volledig gezichtsscherm bij het hanteren van vloeibare stikstof.
7. Breng de snap-frozen monsters over in een vriezer van -80 °C voor langdurige opslag voorafgaand aan de downstream RNA-isolatie en transcriptoomsequencing-pijplijn. Dit is een veilig punt waarop het experiment kan worden gepauzeerd en later opnieuw gestart.
   OPMERKING: De workflow die in dit protocol wordt gepresenteerd, richt zich volledig op het genereren van computationele weerstandshandtekeningen. Er werd geen onafhankelijke experimentele validatie, zoals real-time quantitative PCR (RT-qPCR), uitgevoerd op de belangrijkste differentieel expressieve genen die via deze specifieke pijplijn waren geïdentificeerd.

2. RNA kwaliteitscontrole en bibliotheekvoorbereiding

1. Beoordeel de integriteit van het RNA met behulp van een microfluïdisch capillairelektroforesesysteem. Voor deze representatieve workflow worden RNA-monsters opgenomen met een RNA-integriteitsgetal (RIN) ≥ 6,0, een A260/280-verhouding tussen 1,8 en 2,1, en geen zichtbare degradatiepiek. Noteer de gemeten RIN- en zuiverheidsverhoudingen voor elk monster voordat je de bibliotheekvoorbereiding maakt.
2. Voer 1 μg totaal RNA per monster in voor de bibliotheekvoorbereiding. Zuiver het mRNA uit het totale RNA met behulp van poly-T oligo-gehechte magnetische kralen om te verrijken voor polyA-staartige transcripten.
3. Fragmenteer het verrijkte mRNA met behulp van divalente kationen. Incubeer het mengsel bij 94 °C gedurende 15 minuten in een 5X eerste-streng synthese-reactiebuffer.
4. Synthetiseer het eerste-streng cDNA met behulp van willekeurige hexamerprimers en een omgekeerde transcriptase zonder RNase H-activiteit.
5. Degradeer de RNA-templatestreng met behulp van RNase H. Synthesiseer het tweede streng cDNA met behulp van DNA-polymerase I en dNTP's in een 20 μL reactiesysteem.
6. Incubeer de tweede-streng synthesereactie bij 16 °C gedurende 1 uur. Centrifugeer het reactiemengsel kort op 2.000 x g om de vloeistof op de bodem van de buis op te vangen.
7. Zet de overige overhangen om in stompe uiteinden met behulp van exonuclease- en polymerase-activiteiten. Adenylatie van de 3'-uiteinden van de DNA-fragmenten en ligate-adaptors met haarspeldlusstructuren om voor te bereiden op hybridisatie.
8. Zuiver de bibliotheekfragmenten met magnetische vaste-fase reversibele immobilisatiekralen om cDNA-fragmenten van 370–420 bp lang bij voorkeur te selecteren.
9. Voer ethanolwasbeurten uit tijdens de kraalzuivering. Centrifugeer de buizen op 2.000 x g gedurende 30 seconden om eventueel restanten van ethanol te verzamelen en volledig te verwijderen vóór de uiteindelijke elusie.
10. Voer PCR-amplificatie uit met behulp van een high-fidelity DNA-polymerase, universele PCR-primers en monster-specifieke indexprimers.
11. Voer het PCR-thermisch profiel uit met een initiële denaturatie bij 98 °C gedurende 30 seconden. Daarna volgt 12 cycli van 98 °C gedurende 10 seconden, 60 °C voor 30 seconden en 72 °C voor 30 seconden, eindigend met een laatste verlenging bij 72 °C voor 5 minuten.
12. Zuiver de PCR-producten opnieuw met behulp van de magnetische kralen. Pas dezelfde centrifugatieparameters toe uit stap 2.9 om de uiteindelijke bibliotheek te verkrijgen.
13. Kwantificeer de initiële bibliotheekconcentratie met behulp van een fluorometer. Verdunn de uiteindelijke bibliotheek tot een concentratie van 1,5 ng/μL.
14. Meng de verdunde bibliotheek grondig. Centrifugeer het mengsel op 10.000 x g gedurende 1 minuut bij 4 °C om eventueel achtergebleven vuil te verwijderen vóór de eindanalyse.
15. Beoordeel de insertgrootte van de bibliotheek met behulp van het microfluidische capillairelektroforesesysteem.
16. Kwantificeer de effectieve concentratie van de bibliotheek nauwkeurig via real-time kwantitatieve PCR (qRT-PCR) nadat bevestigd is dat de insertgrootte aan de verwachtingen voldoet. Zorg ervoor dat de concentratie hoger is dan 1,5 nM om bibliotheekstabiliteit en sequencingkwaliteit te garanderen.
    OPMERKING: Dit is een veilig punt waarop het experiment kan worden gepauzeerd. De voorbereide bibliotheken kunnen worden opgeslagen bij -20 °C totdat geclustering en sequencing plaatsvinden.

3. Clustering en transcriptoomsequencing

1. Voer de clustering van de indexgecodeerde samples uit op een geautomatiseerd clustergeneratiesysteem. Gebruik een commerciële paired-end cluster kit volgens de instructies van de fabrikant.
2. Sequentie van de bibliotheekvoorbereidingen op een high-throughput sequencingplatform na succesvolle clustergeneratie. Genereer 150 basisparen (bp) paired-end reads.

4. Kwaliteitscontrole van gegevens en leesmapping

1. Beoordeel de kwaliteit van de ruwe data (FASTQ-formaat) met behulp van fastp v0.23.2 voor rauw-read kwaliteitscontrole en filtering. Noteer de commandoregelparameters in een analyselogboek. In deze workflow werden schone reads gegenereerd door adapter-bevattende reads, reads met poly-N sequenties en low-quality reads met identieke filterinstellingen over samples te verwijderen. Een representatief commando met gepaard einde is beschikbaar in Supplementary File 1.
2. Verwerk de ruwe reads via een geautomatiseerde preprocessingsoftware. Verkrijg schone reads door reads met adapters, reads met poly-N sequenties en low-quality reads te verwijderen. Gebruik identieke filterparameters voor alle samples en registreer na het filteren het behouden leesnummer, Q20, Q30 en GC-inhoud.
3. Bereken de Q20-, Q30- en GC-inhoud van de schone data. Definieer potentiële batchvariabelen vóór downstream analyse, waaronder datum van monsterverzameling, datum van RNA-extractie, batch voor bibliotheekvoorbereiding, sequencing lane en sequencing run.
4. Evalueer batch-effecten met PCA en steekproef-tot-monster correlatieanalyse met behulp van genormaliseerde expressiewaarden. Als de monsters voornamelijk clusteren op basis van technische variabelen in plaats van klinische toestand, documenteer dan de getroffen variabele en neem deze op als een covariaat in de differentiële expressieontwerpformule of pas een gevestigde batch-aanpassingsmethode toe vóór downstream visualisatie.
5. Verkrijg het referentiegenoom (Homo sapiens, GRCh38) en de bijbehorende Ensembl-release-bestanden met 109 genannotatiebestanden voor leesuitlijning.
6. Bouw de index van het referentiegenoom met behulp van HISAT2 v2.0.5.
7. Lijn de gepaarde eind schone reads uit op het referentiegenoom met HISAT2 v2.0.5. Gebruik deze splice-bewuste aligneringsmethode om een database van splice-verbindingen te genereren op basis van het annotatiebestand van het genmodel.

5. Nieuwe transcriptvoorspelling en kwantificatie van genexpressie

1. Stel de toegewezen reads van elk sample samen met StringTie v1.3.3b in een referentiegebaseerde aanpak. Gebruik dit hulpmiddel om volledige transcripten te assembleren en te kwantificeren die meerdere splice-varianten voor elke genlocus vertegenwoordigen.
2. Tel het aantal reads dat aan elk gen is gekoppeld met featureCounts v1.5.0-p3. Gebruik de resulterende ruwe gehele leestellingmatrix als input voor downstream differentiële expressieanalyse.
3. Configureer featureCounts v1.5.0-p3 voor paired-end sequencingdata met de paired-end optie (bijv. -p). Verskaf het gedownloade GRCh38 GTF-annotatiebestand om de juiste genomische kenmerken grenzen te definiëren.
4. Bereken de fragmenten per kilobase van transcript per miljoen gemapte reads (FPKM) voor elk gen. Gebruik FPKM-waarden alleen voor beschrijvende visualisatie, PCA, warmtekaartweergave en verkennende expressiesamenvattingen; gebruik geen FPKM-waarden als invoermatrix voor DESeq2 differentiële expressietests.

6. Analyse van differentiële genexpressie

1. Voer differentiële expressieanalyse uit tussen de nieuw gediagnosticeerde en teruggevallen groepen met behulp van R v3.5.0 en het DESeq2 R-pakket v1.20.0. Importeer de ruwe leestellingmatrix die in stap 5.2 is gegenereerd in de R-omgeving en bewaar FPKM-waarden alleen voor visualisatie en verkennende analyses.
2. Bouw het gespecialiseerde datasetobject dat vereist is door het analysepakket. Voer het specifieke commando uit (bijv. DESeqDataSetFromMatrix()) om de teldatamatrix te binden met de bijbehorende voorbeeldmetadata-tabel.
3. Definieer de experimentele ontwerpformule binnen het softwareobject. Specificeer de klinische toestand (nieuw gediagnosticeerd versus terugval) als primaire variabele voor vergelijking (bijv. ontwerp = ~ aandoening). Als een technische batchvariabele in stap 4.3 wordt geïdentificeerd en niet volledig verward is met de klinische toestand, neem deze dan op in de ontwerpformule (bijv. ontwerp = ~ batch + conditie).
4. Voer de kernfunctie van differentiaalexpressieanalyse uit (bijv. DESeq()). Laat de software automatisch schatting van de groottefactor, spreidingschatting en negatieve binomiale Wald-test passen²² uitvoeren.
5. Haal de resultaattabel uit met behulp van de resultaatextractiefunctie (bijv. resultaten()). Specificeer het contrastargument om de exacte vergelijking te definiëren (teruggevallen versus nieuw gediagnosticeerd).
6. Pas de resulterende P-waarden aan om de valse ontdekkingsfrequentie te controleren. Maak gebruik van de geïntegreerde Benjamini- en Hochberg-procedure die automatisch wordt toegepast door het softwarepakket²³.
7. Filter de tabel met geëxtraheerde resultaten om de significant differentieel tot expressie gebrachte genen (DEGs) te isoleren. Ken elk gen met een aangepaste P-waarde < 0,05 en een absolute logaritmische_2-voudige verandering > 1 toe als significant differentieel uitgedrukt.

7. Genontologie (GO) Verrijkingsanalyse

1. Voer Gene Ontology (GO) verrijkingsanalyse uit van de geïdentificeerde DEG's met behulp van clusterProfiler v3.8.1 en org. Hs.eg.db v3.6.0. Voer de lijst van Entrez-gen-ID's in die overeenkomen met de significante DEG's die in stap 6.7 zijn geïdentificeerd.
2. Voer de GO-verrijkingsfunctie uit (bijv. enrichGO()). Specificeer de vereiste parameters, inclusief de juiste achtergrondorganismendatabase (bijv. OrgDb = org. Hs.eg.db), het specifieke ontologiedomein (Biologisch Proces, Cellulair Component of Moleculaire Functie) en de aangepaste P-waarde cutoff (0,05).
3. Zorg ervoor dat het algoritme de benodigde correcties toepast tijdens de verrijkingsberekening. Bevestig dat de software intern corrigeert voor genlengtebias en de P-waarden aanpast met behulp van de Benjamini- en Hochberg-methode²⁴.
4. Beschouw GO-termen met een gecorrigeerde P-waarde kleiner dan 0,05 als aanzienlijk verrijkt. Genereer een dotplot of staafdiagram met behulp van de geïntegreerde visualisatiefuncties van het pakket om de beste verrijkte GO-termen weer te geven.

8. Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Enrichment Analysis

1. Maak gebruik van een uitgebreide databasebron die gewijd is aan het begrijpen van functies van het biologische systeem op hoog niveau om de ontregelde paden te identificeren. Maak dezelfde lijst van belangrijke DEG Entrez-ID's als gebruikt in stap 7.1.
2. Voer de KEGG-verrijkingsfunctie uit (bijv. enrichKEGG()) binnen het functionele annotatiepakket R.
3. Definieer de kritieke parameters binnen de functieaanroep. Stel de organismecode strikt in op mens (bijv. organisme = 'heeft') en definieer de P-waarde aanpassingsmethode (bijv. pAdjustMethod = 'BH').
4. Extraheer de statistisch significante KEGG-routes. Filter de uitvoer zodat alleen die paden behouden blijven die een gecorrigeerde P-waarde kleiner dan 0,05 aantonen.
5. Visualiseer de meest verrijkte KEGG-paden. Gebruik de geïntegreerde plotfuncties (bijv. dotplot()) om de statistische significantie en gentellingen die bij elk pad horen in kaart te brengen.

9. Genensetverrijkingsanalyse (GSEA)

1. Bereid de vooraf gerangschikte genenlijst voor die nodig is voor de analyse. Bereken de rangschikkingsmetriek voor alle tot expressie gebrachte genen met behulp van de getekende -log10(P-waarde) vermenigvuldigd met het teken van de log 2-voudige verandering afgeleid uit de differentiële expressieanalyse.
2. Start een lokale installatie van de Broad Institute GSEA-software v4.2.3. Voer de nieuw gegenereerde, vooraf gerangschikte genenlijst in in de software-interface²⁵.
3. Download de vereiste vooraf gedefinieerde gensets. Verkrijg de Gene Ontology (GO) en Kyoto Encyclopedia of Genes and Genomes (KEGG) datasets uit de Molecular Signatures Database (MSigDB, versie 7.5.1)²⁶.
4. Configureer de softwareparameters om de statistische verrijkingstest uit te voeren. Stel het aantal permutaties in op 1.000 en selecteer het permutatietype als 'gene_set'.
5. Voer het analyse-algoritme uit om te bepalen of de vooraf gedefinieerde gensets een statistisch significant, overeenkomend verschil laten zien tussen de nieuw gediagnosticeerde en terugvallende biologische toestanden.
6. Evalueer de statistische significantie van de gegenereerde verrijkingsprofielen. Definieer significante gensets met strikte drempels: een genormaliseerde verrijkingsscore (NES) absolute waarde > 1,0, een nominale P-waarde < 0,05 en een valse ontdekkingsrate (FDR) q-waarde < 0,25.

10. Eiwit-eiwitinteractie (PPI) netwerkanalyse

1. Toegang tot de STRING-database voor bekende en voorspelde eiwit-eiwitinteracties. In deze workflow werd de PPI-analyse uitgevoerd met STRING v11.5²⁷.
2. Voer de lijst met Entrez-gen-ID's of officiële gensymbolen in voor de significant differentieel tot expressie gebrachte genen (geïdentificeerd in stap 6.7) in de database-zoekinterface. Selecteer Homo sapiens als doelorganisme.
3. Configureer de netwerkconstructieparameters om te zorgen dat hoogwaardige interacties worden opgehaald. Stel de minimaal vereiste interactiescore in op een hoge betrouwbaarheidsdrempel (score > 0,700).
4. Exporteer de resulterende interactienetwerkgegevens naar een lokale map. Sla de interactiekaart op als een standaard tabelbestand (bijvoorbeeld TSV-formaat).
5. Importeer de geëxporteerde interactiegegevens in Cytoscape v3.9.1 voor netwerkvisualisatie en analyse²⁸.
6. Filter het gebouwde netwerk om de visualisatiehelderheid te verbeteren en belangrijke regelgevende hubs te markeren. Verwijder alle losgekoppelde knooppunten of wees-genen die geen continue interacties vertonen die aan de vastgestelde betrouwbaarheidsdrempel voldoen.

Resultaten

Validatie van klinische cohorten en sequencing

De succesvolle uitvoering van het upstream RNA-extractie- en bibliotheekvoorbereidingsprotocol (Figuur 1) werd bevestigd door sequencing yield- en kwaliteitsmetrics. In deze representatieve dataset leverden beenmergmonsters van vijf nieuw gediagnosticeerde AML-patiënten en vier teruggevallen AML-patiënten gemiddeld ongeveer 6,0 GB ruwe data per monster op. De kwaliteitscontrole (Tabel 2) bevestigde dat de basiskwaliteit en leesdiepte voldoen aan de drempels die vereist zijn voor downstream bio-informatische analyse9. Lage RNA-integriteit (bijvoorbeeld RIN < 6,0), lage mappingsnelheden of hoge transcriptdegradatiebias zouden suboptimale inputkwaliteit vertegenwoordigen en de betrouwbaarheid van downstream differentiële expressieanalyse aantasten.

Globale transcriptomische variantie en PCA

Om de globale transcriptomvariantie te beoordelen en de klinische groepering te inspecteren, werd PCA uitgevoerd op de genormaliseerde expressiegegevens. In deze representatieve dataset vertoonden de nieuw gediagnosticeerde en terugvallende groepen een scheiding in tweedimensionale ruimte (Figuur 2A)²⁰, waarbij PC1 en PC2 respectievelijk 23,82% en 18,75% van de totale variantie uitmaakten. De Venn-diagrammen in Figuur 2B,C bieden een aanvullende beschrijvende samenvatting van genen die zijn gedetecteerd bij monsters binnen de nieuw gediagnosticeerde en terugvallende groepen, ter ondersteuning van reproductieblijkheidscontroles op monsterniveau vóór downstream differentiële expressieanalyse. Omdat de cohort klein en niet gepaard was, werd PCA-scheiding geïnterpreteerd als een illustratieve workflow-output in plaats van als definitief bewijs van ziektetoestand-specifieke biologie.

Differentiële expressiegen (DEG) Analyse

Door de vastgestelde protocoldrempels (|log₂FC| ≥ 1 en aangepaste P-waarde ≤ 0,05) toe te passen op de DESeq2-output werden 2.025 DEGs vastgesteld, bestaande uit 772 upregulated en 1.253 downregulated genen in de terugkerende groep (Figuur 3A). Kandidaattranscripties met grote variatie omvatten FOXC1 (log₂FC = 7,55, P = 4,92 x ^10-5), HOXA11 (log₂FC = 7,76), HOXA11-AS (log₂FC = 7,23) en AXL (log₂FC = 3,50), naast downregulated RHOB, PTX3 en CXCL8. Bestaande literatuur koppelt verschillende van deze genen aan AML-stamvorm, signalering of therapierespons^13,29; de huidige workflow identificeert ze echter alleen als terugvalgerelateerde kandidaat-transcripties. Elke definitieve mechanistische rol in klinische resistentie vereist daaropvolgende onafhankelijke functionele validatie.

Functionele en padverrijking (GO, KEGG en GSEA)

Het functionele annotatieprotocol koppelde de DEG's aan bredere biologische systemen. GO-analyse identificeerde verrijking van termen gerelateerd aan kleine GTPase-gemedieerde signaaltransductie, metaaliontransport en chromatineassemblage (Figuur 4A–C). KEGG-routekaartlegging identificeerde associaties met ECM-receptorinteracties en cytokine-cytokine receptorinteracties (Figuur 4D). GSEA toonde verrijking van RNA-biosynthetische processen in de terugvallende groep en verrijking van energiemetabolismeroutes in de nieuw gediagnosticeerde groep (Figuur 5A). Deze verrijkingsresultaten bieden een beschrijvende routekaart van veranderde gensets en moeten worden geïnterpreteerd als hypothese-genererende associaties in plaats van bewezen oorzaken van terugval.

Constructie van het eiwit-eiwitinteractienetwerk (PPI)

Het oorspronkelijke STRING-netwerk bevatte 56 knooppunten en 193 interacties. Na het verwijderen van losgekoppelde of wees-knopen bevatte het weergegeven Cytoscape-subnetwerk 42 knopen en 136 interacties (Figuur 5B). Netwerkmodulaire analyse gaf prioriteit aan TP53, CCL2, CXCL8 en IL6 als centrale wiskundige hubs met het hoogste aantal interacties. Omdat het PPI-netwerk afhankelijk is van database-voorspelde interactiescores (bijv. ATF3-score: 0,982), moet hubidentificatie worden geïnterpreteerd als doelprioritering voor toekomstige empirische studies in plaats van direct bewijs van p53-gemedieerde apoptoseontwijking of andere weerstandsmechanismen.

De ruwe RNA-sequencinggegevens die in dit protocol worden gegenereerd, zijn gedeponeerd in de Figshare-repository en zijn openbaar toegankelijk via de volgende DOI: https://doi.org/10.6084/m9.figshare.30655814. De verwerkte gegevens en bijbehorende analysebestanden zijn opgenomen in het artikel en het aanvullende materiaal. Representatieve commandoregelparameters en analyse-instellingen die worden gebruikt om de computationele workflow te reproduceren, worden geleverd als Supplementary File 1. Alle gegevens die de bevindingen van deze studie ondersteunen, zijn zonder beperkingen beschikbaar.

Patiëntidentificatie	Leeftijd (jaren)	Seks	Moleculaire mutaties	Overleven/Follow-up (maanden)	Klinische status
R_AML_1	70	Mannelijk	FLT3-ITD (+)	22	Overleden
R_AML_2	29	Vrouwelijk	NPM1 (+)	11	Levend
R_AML_3	40	Mannelijk	CEBPA (+)	17	Levend
R_AML_4	55	Vrouwelijk	Drievoudig Negatief*	24	Overleden

Tabel 1: Demografische en klinische kenmerken van patiënten in de terugvallende AML (R_AML) groep. Tabel 1 vat de demografische en klinische kenmerken samen van de terugvallende AML-cohort die in de representatieve analyse werd gebruikt, inclusief klinische kenmerken op patiëntniveau die relevant zijn voor de interpretatie van de transcriptomische workflow.

Voorbeeld	Bibliotheek	Raw_reads	Raw_bases	Clean_reads	Clean_bases	Error_rate	Q20	Q30	GC_pct
AML_1	FRAS25 0244891-1r	48705066	7.31G	47807532	7,17G	0.01	99.35	97.48	47.48
AML_2	FRAS25 0244896-1r	42969940	6,45G	42237962	6,34G	0.01	99.35	97.44	46.74
AML_3	FRAS2502 44906-1r	48738386	7.31G	47744462	7,16G	0.01	99.36	97.48	47.28
AML_4	FRAS250 244915-1r	48723650	7.31G	47688240	7,15G	0.01	99.29	97.26	47.45
AML_5	FRAS2502 44920-1r	49508198	7,43G	47740308	7,16G	0.01	99.37	97.53	47.73
R_AML_1	FRAS2502 44892-1r	47879408	7,18G	46671584	7.0G	0.01	99.39	97.49	47.63
R_AML_2	FRAS2502 70005-1r	47657378	7,15G	46957882	7.04G	0.01	99.39	97.49	50.5
R_AML_3	FRAS250 405722-1r	58754766	8,81G	56867112	8,53G	0.01	99.38	97.42	46.52
R_AML_4	FRAS2502 44902-1r	48491122	7,27G	47469334	7.12G	0.01	99.23	97.21	46.43

Tabel 2: Samenvatting van de gegevenskwaliteit. Tabel 2 geeft sequencing-kwaliteitsmetrics voor elk monster weer, inclusief leesrendement, basiskwaliteit, GC-inhoud en mapping-gerelateerde kwaliteitscontrole-informatie die wordt gebruikt om te bepalen of monsters geschikt zijn voor downstream-analyse.

Figuur 1: Workflow van het protocol. De workflow vat de belangrijkste experimentele en computationele fasen samen, waaronder klinische monsterverzameling, RNA-kwaliteitscontrole, bibliotheekvoorbereiding en -sequencing, leesverwerking en -uitlijning, transcriptkwantificatie, differentiële expressieanalyse, GO/KEGG-verrijking, GSEA en PPI-netwerkconstructie. Klik hier om een grotere versie van deze figuur te bekijken.

Figuur 2: Kwantitatieve analyse van steekproeven. (A) Principal component analysis (PCA) werd uitgevoerd om intergroepverschillen en reproduceerbaarheid binnen de groep te evalueren. PCA werd uitgevoerd met lineaire algebraïsche methoden gebaseerd op genormaliseerde genexpressiewaarden over alle monsters. (B, C) Venn-diagrammen die genen tonen die worden gedetecteerd in monsters in respectievelijk de AML- en R_AML groepen. Monster-beperkte regio's geven genen aan die in individuele monsters zijn gedetecteerd, terwijl overlappende gebieden genen vertegenwoordigen die vaak worden gedetecteerd over twee of meer monsters. Klik hier om een grotere versie van deze figuur te bekijken.

Figuur 3: Analyse van differentiële genexpressie. (A) Balkgrafiek die het aantal differentieel tot expressie gebrachte genen (DEGs) tussen vergelijkingsgroepen toont, geïdentificeerd door DESeq2 met drempels van aangepaste P-waarde ≤ 0,05 en |log₂FoldChange| ≥ 1. (B) Vulkaanplot van DEGs. De x-as vertegenwoordigt log₂FoldChange-waarden, en de y-as staat voor -log10(P-waarde). Blauwe stippellijnen geven de drempellijnen aan die worden gebruikt voor DEG-selectie. (C) Hiërarchische clustering-heatmap van LEG's. De x-as geeft voorbeeldnamen aan, en de y-as toont genormaliseerde expressiewaarden van de LEG's. Klik hier om een grotere versie van deze figuur te bekijken.

Figuur 4: Analyse van functioneel verrijking van differentieel tot expressie gebrachte genen. (A) GO verrijkingsbalk-plot. De x-as stelt GO-termen weer en de y-as toont verrijkingssignificantie, uitgedrukt als -log10(padj). Kleuren staan voor BP (Biologisch Proces), CC (Cellulaire Component) en MF (Moleculaire Functie). (B) GA verrijkingsbubbelplot. De x-as geeft de verhouding van DEG's aan die aan elke GO-term zijn geannoteerd ten opzichte van het totale aantal DEG's, en de y-as geeft de GO-termen aan. De bubbelgrootte komt overeen met het aantal geannoteerde genen, en kleurgradiënten geven de verrijkingssignificantie weer aan. (C) KEGG verrijkingsbalkplot. De x-as vertegenwoordigt KEGG-paden, en de y-as geeft de verrijkingsbetekenis aan. (D) KEGG verrijkingsbubbelplot. De bubbelgrootte geeft het aantal geannoteerde genen aan, en kleurgradiënten weerspiegelen de verrijkingsbetekenis. Klik hier om een grotere versie van deze figuur te bekijken.

Figuur 5: GSEA-verrijking en eiwit-eiwitinteractie (PPI) netwerkanalyse. (A) Bargrafiek die genormaliseerde verrijkingsscores (NES) toont voor geselecteerde significante gensets. Positieve NES-waarden duiden op verrijking in de R_AML-groep, terwijl negatieve NES-waarden verrijking aangeven in de nieuw gediagnosticeerde AML-groep. (B) Eiwit-eiwitinteractie (PPI) netwerk. Elke knoop vertegenwoordigt een eiwit, en elke rand duidt op een interactie tussen verbonden eiwitten. Klik hier om een grotere versie van deze figuur te bekijken.

Discussie

Kritieke stappen in het protocol

De succesvolle uitvoering van deze bio-informatische workflow is afhankelijk van verschillende kritieke stappen. Ten eerste zijn het onmiddellijke bevriezen en de juiste lysatie van het beenmergaspirat (stap 1.6) van groot belang, aangezien de micro-omgeving van het beenmerg rijk is aan ribonucleasen die transcriptomische integriteit^{snel kunnen} aantasten. Tijdens de computationele fase is de juiste configuratie van de experimentele ontwerpformule binnen het DESeq2-pakket (stap 6.3) cruciaal voor een nauwkeurige differentiële expressie, vooral wanneer de klinische toestand (nieuw gediagnosticeerd versus recidiverend) wordt gecontrasteerd terwijl potentiële verstorende variabelen worden meegenomen. Ten slotte is het toepassen van rigoureuze false discovery rate (FDR) drempels tijdens Gene Set Enrichment Analysis (GSEA) (stap 9.6) een cruciaal statistisch controlepunt om overinterpretatie van vals-positieve functionele netwerken te voorkomen.

Aanpassingen en probleemoplossing

Een veelvoorkomende uitdaging bij deze methode is de aanwezigheid van batch-effecten, die vaak optreden wanneer klinische monsters worden verzameld en gesequenced over langere tijdsperioden. Batchvariabelen moeten vóór de analyse worden gedefinieerd, waaronder de datum van monsterverzameling, de datum van RNA-extractie, de batch van bibliotheekvoorbereiding, de sequencing-lane en de sequencing-run. Als PCA of steekproefcorrelatieanalyse clustering aantoont op basis van sequencingdatum of een andere technische variabele in plaats van klinisch fenotype, dienen gebruikers het protocol aan te passen door de batchvariabele op te nemen in de differentiële expressieontwerpformule wanneer statistisch haalbaar of door batchcorrectie-algoritmen toe te passen, zoals ComBat of SVA, vóór visualisatie³¹. Bij toepassing van dit protocol op volbloedaspiraten in plaats van beenmergaspiraten, is een essentiële aanpassing de opname van een globine-mRNA-depletiestap tijdens de bibliotheekvoorbereiding om te voorkomen dat zeer overvloedige globinetranscripten de leesdiepte van sequencing monopoliseren. Softwareversies en hoofdparameters voor de representatieve workflow werden als volgt aangevuld: 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, gepaarde 150 bp sequencing, GSEA v4.2.3 met 1.000 gensetpermutaties, MSigDB v7.5.1, STRING v11.5 met interacties met hoge betrouwbaarheid, en Cytoscape v3.9.1. Representatieve commandoregelparameters en analyse-instellingen zijn beschikbaar in Supplementair Bestand 1.

Beperkingen van de methode

Hoewel het protocol volledig is, kent het inherente methodologische beperkingen. Ten eerste gebruikt het bulk RNA-sequencing, dat het gemiddelde transcriptomische profiel van de gehele beenmergaspirate vastlegt en geen ruimtelijke resolutie heeft voor één cel. Daarom kan de workflow niet bepalen of een opgereguleerde terugval-geassocieerde signatuur afkomstig is van leukemie-stamcellen, stromale cellen, immuuncellen of veranderingen in celtype-samenstelling³². Ten tweede is de representatieve dataset klein (n = 9) en ongepaard, wat de statistische robuustheid beperkt en definitieve causale inferentie voorkomt. Ten derde is de workflow volledig in silico. Het genereert kandidaat-regulatorische hubs en signaalroutes, maar kan hun functionele noodzaak bij chemoresistentie niet onafhankelijk valideren zonder orthogonale in vitro of in vivo experimentele validatie.

Recente enkelcel- en enkelcel-genomische studies hebben het AML-referentiekader uitgebreid door cel-toestandheterogeniteit, klonale architectuur en therapie-geassocieerde evolutie op te lossen bij hogere resolutie 33,34,35,36. Deze benaderingen zijn complementair aan de hier beschreven bulk RNA-seq workflow: bulk sequencing biedt een praktische en kosteneffectieve screeningsstrategie voor transcriptomische signaturen op cohortniveau, terwijl single-cell en multi-omic methoden in vervolgstudies kunnen worden gebruikt om kandidaatsignalen toe te wijzen aan specifieke kwaadaardige of micro-omgevingscelpopulaties.

Betekenis met betrekking tot bestaande methoden

Ondanks deze beperkingen biedt deze transcriptomische pijplijn voordelen ten opzichte van alternatieve diagnostische en analytische technieken. Traditionele klinische beoordelingen van AML-terugval zijn vaak afhankelijk van gerichte multiplex qPCR-panelen of standaard flowcytometrie. Hoewel nuttig voor snelle diagnostiek, zijn deze gerichte methoden beperkt door vooraf gedefinieerde probes en kunnen ze alleen bekende weerstandsmarkers¹⁹ evalueren. Door gebruik te maken van onbevooroordeelde, genoombrede transcriptomsequencing in combinatie met netwerkanalyse, kan dit protocol nieuwe transcripten en systeembrede associaties nomineren die bestaande gerichte methoden mogelijk over het hoofd zien.

Belang en Potentiële Toepassingen

De methodologie die in dit protocol wordt beschreven, is relevant voor translationele hematologie en gepersonaliseerde geneeskunde omdat deze transcriptomische signaturen die met terugval geassocieerd zijn kan prioriteren voor aanvullend onderzoek. Een mogelijke downstream toepassing is het aanwijzen van oppervlakte-antigenen of immuunontwijkingsroutes die ontstaan tijdens een terugval. Dergelijke kandidaten kunnen het ontwerp van toekomstige validatiestudies informeren en, indien experimenteel bevestigd, kunnen bijdragen aan de ontwikkeling van de volgende generatie immunotherapieën, waaronder CAR-T of CAR-NK celstrategieën³⁷. Deze translationele toepassingen moeten worden beschouwd als hypothesegenererend in plaats van als vastgestelde conclusies uit de huidige dataset.

Materialen

Lijst van materialen gebruikt in dit artikel

Naam	Bedrijf	Catalogusnummer	Opmerkingen
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 or institution-specific system ID
clusterProfiler (Software)	Bioconductor	RRID:SCR_016884	v3.8.1
Cytoscape (Software)	Cytoscape Consortium	RRID:SCR_003032	v3.9.1
DESeq2 (Software)	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 (Software)	Bioconductor	RRID:SCR_012802	v3.22.5
fastp (Software)	OpenGene	RRID:SCR_016962	v0.23.2
featureCounts / Subread (Software)	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 (Software)	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 or laboratory-specific kit
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; service-provider instrument ID
org.Hs.eg.db (Annotation package)	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	-	Project No. X101SC25054246-Z01-J003
STRING database	STRING Consortium	RRID:SCR_005223	v11.5
StringTie (Software)	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