Transkriptomisches Profilieren und bioinformatische Analyse von Knochenmarkproben zur Identifizierung von Chemotherapieresistenzsignaturen bei akuter myeloider Leukämie

Zusammenfassung

Dieses Protokoll bietet einen standardisierten bioinformatischen Workflow zur Analyse transkriptomischer Veränderungen bei akuter myeloischer Leukämie (AML). Ziel ist es, neu diagnostizierte und rückfällige Knochenmarkproben zu vergleichen und molekulare Signaturen, die mit Chemotherapieresistenz und dem Krankheitsverlauf verbunden sind, für weitere Untersuchungen zu priorisieren.

Zusammenfassung

Die akute myeloide Leukämie (AML) ist eine hoch heterogene hämatologische Malignität, bei der Rückfälle und erworbene Chemoresistenz Hauptursachen für das Therapieversagen bleiben. Dieser Artikel stellt ein bioinformatisches Protokoll zur transkriptomischen Analyse von Knochenmarkaspirationen vor. Das Hauptziel des Protokolls ist es, einen standardisierten Arbeitsablauf bereitzustellen, um molekulare Signaturen zu identifizieren, die mit dem Krankheitsverlauf und der Therapieresistenz bei rezidieren AML verbunden sind. Die Pipeline beschreibt die rechnerischen Verfahren zum Vergleich ungepaarter Knochenmarkproben, demonstriert durch Sequenzierungsdaten von fünf neu diagnostizierten Fällen und vier rückfälligen Fällen. Diese Methode beschreibt die wesentlichen Schritte zur Verarbeitung von RNA-Sequenzierungsdaten, zur Durchführung differenzieller Genexpressionsanalysen und zur Durchführung von nachgelagerten funktionellen Evaluationen. Durch die Anwendung dieses Workflows wurden in diesem repräsentativen Datensatz 2.025 differenziell exprimierte Gene (DEGs), darunter FOXC1, HOXA11, HOXA11-AS und AXL, als Kandidatentranskripte identifiziert, die mit einem Rückfall assoziiert sind. Funktionelle und Netzwerkanalysen priorisierten zudem Gensätze und Interaktionshubs im Zusammenhang mit der Signalübertragung kleiner GTPase, entzündlicher Signalübertragung, extrazellulären Matrix-Interaktionen und RNA-Biosynthetischen Prozessen. Insgesamt bietet diese Methodik eine reproduzierbare Rechenpipeline zur Abbildung transkriptomischer Signaturen, die mit rückfälligem AML assoziiert sind, und zur Generierung von Hypothesen, die eine anschließende experimentelle Validierung erfordern.

Einleitung

Akute myeloide Leukämie (AML) ist eine Gruppe klonaler bösartiger Neoplasien, die aus hämatopoietischen Stamm-/Vorläuferzellen stammen und durch eine abnormale Proliferation unreifer myeloider Zellen im Knochenmark und Unterdrückung der hämatopoetischen Differenzierung^{gekennzeichnet sind.} Obwohl die derzeitige Standard-Induktionsjemotherapie (wie Cytarabin in Kombination mit Anthracyklinen) bei den meisten Patienten eine vollständige Remission auslösen kann, bleibt die Rückfallrate bei bis zu 50 % bis 70 % hoch, und die Prognose für rückfällige Patienten ist signifikant beeinträchtigt ^4,5,6. Die erworbene Chemotherapieresistenz ist mit dem Versagen der AML-Behandlung assoziiert, was eine eingehende Analyse der mit diesem Prozess verbundenen molekularen Signaturen erforderlich macht, um Therapiestrategien und Patientenüberlebensraten zu verbessern ^7,8.

In der breiteren Literatur zeigen bestehende Studien, dass die Chemotherapieresistenz bei AML nicht auf die Hochregulation von Medikamenten-Effluxpumpen oder einen abnormalen Stoffwechsel beschränkt ist, sondern auch mit der Überlebensaufrechterhaltung von Leukämie-Stammzellen (LSC), der Bildung von epithelial-mesenchymalen Übergangsformen (EMT)-ähnlichen Phänotypen innerhalb der hämatologischen Nische und Interaktionen mit der Knochenmarkmikroumgebung^{9, 10,11}. So zeigen beispielsweise LSC-Populationen eine hohe Selbsterneuerungskapazität und Ruhe, was mit einer inhärenten Resistenz gegen zellzyklusspezifische Chemotherapeutika einhergeht. Zusätzlich wurde die Hochregulierung von Rezeptor-Tyrosinkinasen, wie AXL, mit Resistenzen in FLT3-ITD+AML assoziiert, die parallel zur Aktivierung der PI3K/AKT- und MAPK-Signalwege sowie verbesserten anti-apoptotischen Fähigkeiten^{auftritt 13,14}. Metabolische Umprogrammierung und epigenetische Remodellierung wurden ebenfalls als wichtige regulatorische Achsen bei der Resistenzbildung anerkannt. Belege deuten darauf hin, dass AML-Zellen während eines Rückfalls durch Chemotherapie verursachten oxidativen Stress und DNA-Schäden durch verstärkte oxidative Phosphorylierungsaktivität (OXPHOS), modulierte NAD⁺/NADH-Verhältnisse und veränderte Histonmodifikationszustände 15,16,17 anpassen können. Entzündungsfaktoren im Knochenmarkmikroumfeld wie IL-6 und CXCL8 sind ebenfalls mit LSC-Überleben und Chemotherapieresistenz assoziiert und treten häufig in Koordination mit der Aktivierung der STAT3/NF-κB-Signalwege^11,18 auf.

Trotz dieser anerkannten Mechanismen bleiben die transkriptomischen Veränderungen, die mit AML-Rückfall und Chemoresistenz assoziiert sind, unvollständig charakterisiert, insbesondere wenn metabolische, epigenetische und knochenmarkbezogene Mikroumgebungssignale gemeinsam in klinischen Proben bewertet werden. Dieser Workflow adressiert diesen Bedarf, indem er Kandidaten-DEGs, Signalwege und regulatorische Netzwerke priorisiert, die mit dem Übergang von der Erstdiagnose zum klinischen Rückfall verbunden sind. Der Ansatz integriert differentielle Expressionsanalyse, Genset-Anreicherungsanalyse (GSEA) und den Aufbau von Protein-Protein-Interaktionsnetzwerken (PPI), um systemweite transkriptionelle Neuprogrammierung abzubilden und Kandidaten für eine anschließende mechanistische Validierung zu nominieren.

Das übergeordnete Ziel dieser Methode ist es, eine standardisierte, reproduzierbare bioinformatische Pipeline zum Vergleich der Bulk-Transkriptome neu diagnostizieller mit rezidivierten AML-Knochenmarkproben zu präsentieren. Die Begründung für den Einsatz dieser in-silico-Technik ist ihre Fähigkeit, unvoreingenommene, genomweite transkriptionelle Ereignisse zu erfassen und über die Grenzen von Einzelweganalysen hinauszugehen und komplexe, multidimensionale regulatorische Netzwerke zu priorisieren. Diese Technik bietet erhebliche Vorteile gegenüber alternativen Methoden wie Mikroarrays oder gezielten multiplexen qPCR-Panels, da sie einen höheren dynamischen Bereich, die Fähigkeit zur Erkennung neuartiger Transkripte und eine präzise Quantifizierung der Genexpression ohne die Einschränkungen vorgefertigter Sonden^bietet. Um festzustellen, ob diese Methode für ihre Anwendung geeignet ist, sollten Leser beachten, dass diese Pipeline speziell für Forscher konzipiert ist, die Massen-RNA-Sequenzierungsdaten von gepaarten oder ungepaarten klinischen Kohorten, wie Gewebeaspirationen, verarbeiten. Es eignet sich zur Identifizierung breit widerstandsassoziierter Signaturen und Kandidaten-Regulationsnetzwerke, während Forscher, die zelltypspezifische oder räumliche Auflösung benötigen, komplementäre Einzelzell- oder räumliche Sequenzierungs-Workflows verwenden müssten. Letztlich ermöglicht dieses rechnergestützte Protokoll die Priorisierung von Kandidatengenen, Signalwegen und regulatorischen Netzwerken für spätere experimentelle Untersuchungen.

Protokoll

Alle Methoden zur Verwendung menschlicher Gewebeproben wurden in Übereinstimmung mit institutionellen Richtlinien und der Helsinkier Erklärung (überarbeitet 2013) durchgeführt. Die klinischen Knochenmarkproben wurden mit Zustimmung des institutionellen Ethikausschusses gewonnen (Genehmigungsnummern TY-ZKY2024-116-01 und TY-ZKY2024-116-02).

1. Klinische Probenentnahme und Patientenklassifikation

1. Wählen Sie Knochenmarkaspirationsproben von Patienten aus, bei denen anhand der Klassifikationskriterien der Weltgesundheitsorganisation (WHO) formell akute myeloische Leukämie (AML) diagnostiziert wurden.
2. Anwendung spezifischer Ein- und Ausschlusskriterien bei der Patientenauswahl, um Kohortenhomogenität und Reproduzierbarkeit sicherzustellen. Erwachsene Patienten mit primärer AML einbeziehen und Patienten mit sekundärer AML, akuter promyelocytärer Leukämie oder einer Vorgeschichte anderer Malignitäten ausschließen (Tabelle 1).
3. Weisen Sie die gesammelte Probe der neu diagnostizierten Gruppe zu, wenn der Patient zum Zeitpunkt der ersten klinischen Diagnose mit zuvor unbehandeltem AML vorliegt.
4. Die gesammelte Probe wird der rückfälligen Gruppe zugeordnet, wenn der Patient nach einer dokumentierten vollständigen Remission ein erneutes Auftreten leukämischer Blasten im peripheren Blut oder mehr als 5 % im Knochenmark nachweist.
5. Sammeln Sie die deidentifizierten Restproben des Knochenmarks unmittelbar nach dem routinemäßigen klinischen Knochenmarkaspirationsverfahren.
   HINWEIS: Bei der Demonstration dieses spezifischen Protokolls wurden von September 2024 bis September 2025 neun aufeinanderfolgende Proben (fünf neu diagnostizierte und vier rückfällige) entnommen. Da dabei nur deidentifizierte Restproben verwendet wurden, verzichtete der Ethikausschuss auf die Verpflichtung einer schriftlichen informierten Einwilligung.
6. Verarbeiten Sie das gesammelte Knochenmarkaspirat sofort zur RNA-Erhaltung mit der Standard-Phenol-Guanidinium-Lysemethode²¹.
   1. Übertragen Sie das frisch entnommene Knochenmark-Aspiratium in ein Sammelrohr mit einem Antikoagulans. Schüttle den Schlauch, um das Aspirationsmittel und das Antikoagulans gründlich zu mischen.
   2. Extrahieren Sie ein gemessenes Volumen der antikoagulierten Probe und geben Sie es einem kommerziellen Phenol-Guanidinium-Lysereagenz zu. Halten Sie ein Volumenverhältnis von 3 Teilen Lysereagenz zu 1 Teil Probe ein.
      VORSICHT: Das Phenol-Guanidinium-Lyse-Reagenz enthält giftige und korrosive Chemikalien, die schwere Verbrennungen und Gewebeschäden verursachen können. Führen Sie alle Handhabungen mit Reagenzen in einer chemischen Abzughaube durch, während Sie geeignete persönliche Schutzausrüstung tragen.
   3. Schütteln Sie das Röhrchen kräftig, um die Probe und das Lysereagenz vollständig zu homogenisieren. Stellen Sie sicher, dass die Mischung vollständig vermischt ist, und stellen Sie sicher, dass keine sichtbaren Gerinnsel mehr in der Lösung verbleiben.
   4. Die homogenisierte Probe wird sofort mit einem Snap-Free-Verfahren eingefroren, indem man das Rohr in flüssigen Stickstoff taucht.
      VORSICHT: Flüssigstickstoff ist extrem kalt und kann bei Kontakt starke Erfrierungen verursachen. Tragen Sie kryogene Handschuhe und einen Schutzschild beim Umgang mit flüssigem Stickstoff.
7. Übertragen Sie die Schnappfrorenproben in einen -80 °C Gefrierschrank zur Langzeitlagerung vor der nachgelagerten RNA-Isolierung und der Transkriptom-Sequenzierung. Dies stellt einen sicheren Punkt dar, an dem das Experiment pausiert und später neu gestartet werden kann.
   HINWEIS: Der in diesem Protokoll dargestellte Arbeitsablauf konzentriert sich vollständig auf die Erzeugung von Algorithmenresistenzsignaturen. Eine unabhängige experimentelle Validierung, wie etwa Echtzeit-quantitative PCR (RT-qPCR), wurde an den durch diesen spezifischen Pipeline identifizierten wichtigsten differentiell exprimierten Genen durchgeführt.

2. RNA-Qualitätskontrolle und Bibliotheksvorbereitung

1. Beurteile die RNA-Integrität mit einem mikrofluidischen Kapillarelektrophoresesystem. Für diesen repräsentativen Workflow sollten RNA-Proben mit einer RNA-Integritätszahl (RIN) ≥ 6,0, einem A260/280-Verhältnis zwischen 1,8 und 2,1 sowie ohne sichtbaren Abbau-Peak aufgenommen werden. Erfassen Sie die gemessenen RIN- und Reinheitsverhältnisse für jede Probe vor der Bibliotheksvorbereitung.
2. Geben Sie 1 μg Gesamt-RNA pro Probe für die Bibliotheksvorbereitung ein. Reinigen Sie die mRNA aus der Gesamt-RNA mithilfe von poly-T-oligogebundenen magnetischen Perlen, um für polyA-schwanzige Transkripte zu bereichern.
3. Fragmentieren Sie die angereicherte mRNA mit divalenten Kationen. Inkubieren Sie das Gemisch bei 94 °C für 15 Minuten in einem 5-fachen Erststrang-Synthese-Reaktionspuffer.
4. Synthetisieren Sie die cDNA des ersten Strangs mit zufälligen Hexamerprimern und einer Reverse-Transkriptase ohne RNase-H-Aktivität.
5. Abbau des RNA-Template-Strangs mit RNase H. Synthetisieren Sie die cDNA des zweiten Strangs unter Verwendung von DNA-Polymerase I und dNTPs in einem 20-μL-Reaktionssystem.
6. Inkubieren Sie die Synthesereaktion des zweiten Strangs bei 16 °C für 1 Stunde. Zentrifugiere das Reaktionsgemisch kurz mit 2.000 x g, um die Flüssigkeit am Boden des Rohrs aufzufangen.
7. Wandeln Sie die verbleibenden Überhänge mithilfe von Exonuklease- und Polymeraseaktivitäten in stumpfe Enden um. Adenylieren Sie die 3'-Enden der DNA-Fragmente und Ligatenadapter mit Haarnadelschleifenstrukturen, um sich auf die Hybridisierung vorzubereiten.
8. Reinigen Sie die Bibliotheksfragmente mit magnetischen Festphasen-reversiblen Immobilisierungsperlen, um cDNA-Fragmente von 370–420 bp Länge bevorzugt auszuwählen.
9. Führen Sie während der Perlenreinigung Ethanolwäschen durch. Zentrifugieren Sie die Röhren bei 2.000 x g für 30 Sekunden, um vor der endgültigen Elution den Rest von Ethanol zu sammeln und vollständig zu entfernen.
10. Führen Sie PCR-Amplifikation mit einer hochpräzisen DNA-Polymerase, universellen PCR-Primern und probenspezifischen Indexprimern durch.
11. Führen Sie das PCR-Wärmeprofil mit einer anfänglichen Denaturierung bei 98 °C für 30 Sekunden aus. Es folgt 12 Zyklen mit 98 °C für 10 Sekunden, 60 °C für 30 Sekunden und 72 °C für 30 Sekunden, die mit einer abschließenden Verlängerung bei 72 °C für 5 Minuten enden.
12. Reinigen Sie die PCR-Produkte erneut mit den magnetischen Perlen. Wenden Sie die identischen Zentrifugationsparameter aus Schritt 2.9 an, um die endgültige Bibliothek zu erhalten.
13. Quantifizieren Sie die anfängliche Bibliothekskonzentration mit einem Fluorometer. Verdünnen Sie die Endbibliothek auf eine Konzentration von 1,5 ng/μL.
14. Mische die verdünnte Bibliothek gründlich. Zentrifugieren Sie das Gemisch bei 10.000 x g für 1 Minute bei 4 °C, um eventuelle verbleibende Ablagerungen vor der Endanalyse zu entfernen.
15. Bewerten Sie die Insert-Größe der Bibliothek mit dem mikrofluidischen Kapillarelektrophoresesystem.
16. Quantifizieren Sie die effektive Konzentration der Bibliothek genau mittels Echtzeit-quantitativer PCR (qRT-PCR), nachdem bestätigt wurde, dass die Insert-Größe den Erwartungen entspricht. Stellen Sie sicher, dass die Konzentration höher als 1,5 nM liegt, um die Bibliotheksstabilität und Sequenzierungsqualität zu gewährleisten.
    HINWEIS: Dies stellt einen sicheren Punkt dar, an dem das Experiment pausiert werden kann. Die vorbereiteten Bibliotheken können bei -20 °C bis zur Clusterung und Sequenzierung gespeichert werden.

3. Clustering und Transkriptomsequenzierung

1. Führen Sie das Clustern der indexcodierten Samples in einem automatisierten Cluster-Generierungssystem durch. Verwenden Sie ein kommerzielles Paired-End-Cluster-Set gemäß den Anweisungen des Herstellers.
2. Sequenzieren Sie die Bibliotheksvorbereitungen auf einer Hochdurchsatz-Sequenzierungsplattform nach erfolgreicher Cluster-Generierung. Erzeugen Sie 150 Basepaar-(bp) Paired-End-Reads.

4. Datenqualitätskontrolle und Lesekartierung

1. Bewertung der Qualität der Rohdaten (FASTQ-Format) mit fastp v0.23.2 für die Rohlese-Qualitätskontrolle und Filterung. Zeichnen Sie die Kommandozeilenparameter in einem Analyseprotokoll fest. In diesem Workflow wurden saubere Lesearten erzeugt, indem Adapter-enthaltene Lese, Lesungen mit Poly-N-Sequenzen und minderwertige Lesungen mit identischen Filtereinstellungen über Samples entfernt wurden. Ein repräsentativer Befehl mit gepaartem Ende ist in Supplementary File 1 bereitgestellt.
2. Verarbeiten Sie die Rohlesungen durch eine automatisierte Vorbearbeitungssoftware. Erhalten Sie saubere Lese, indem Sie Lesearten mit Adaptern, Lesearten mit Poly-N-Sequenzen und Lesearten niedriger Qualität entfernen. Verwenden Sie identische Filterparameter für alle Samples und zeichnen Sie die gespeicherten Lesezahlen, Q20, Q30 und GC-Inhalte nach dem Filtern auf.
3. Berechnen Sie den Q20-, Q30- und GC-Inhalt der sauberen Daten. Definieren Sie potenzielle Chargenvariablen vor der nachgelagerten Analyse, einschließlich Probenentnahmedatum, RNA-Extraktionsdatum, Bibliotheksvorbereitungscharge, Sequenzierungsleitung und Sequenzierungslauf.
4. Bewerten Sie Batch-Effekte mittels PCA und Proben-zu-Proben-Korrelationsanalyse mit normalisierten Expressionswerten. Wenn die Stichproben hauptsächlich nach technischen Variablen und nicht nach klinischem Zustand gruppieren, dokumentieren Sie die betroffene Variable und fügen sie als Kovariate in die Differentialexpressions-Designformel ein oder wenden eine etablierte Batch-Anpassungsmethode vor der nachgelagerten Visualisierung an.
5. Erhalte das Referenzgenom (Homo sapiens, GRCh38) und die entsprechenden Ensembl-Release-Dateien mit 109 Genannotationsdateien zur Leseausrichtung.
6. Erstelle den Index des Referenzgenoms mit HISAT2 v2.0.5.
7. Richten Sie die gepaarten End-Clean Reads mit HISAT2 v2.0.5 auf das Referenzgenom aus. Verwenden Sie diesen splice-bewussten Ausrichtungsansatz, um eine Datenbank von Splice-Junctions basierend auf der Genmodell-Annotationsdatei zu erstellen.

5. Neuartige Transkriptvorhersage und Genexpressionsquantifizierung

1. Assembleren Sie die abgebildeten Reads jeder Probe mit StringTie v1.3.3b in einem referenzbasierten Ansatz. Nutzen Sie dieses Werkzeug, um vollständige Transkripte zusammenzustellen und zu quantifizieren, die mehrere Spleißvarianten für jeden Genlokus repräsentieren.
2. Zählen Sie die Anzahl der auf jedes Gen abgebildeten Lesungen mit featureCounts v1.5.0-p3. Verwenden Sie die resultierende rohe ganzzahlige Lesezählmatrix als Eingabe für die nachgeschaltete Differentialausdrucksanalyse.
3. Konfigurieren Sie featureCounts v1.5.0-p3 für Paired-End-Sequenzierungsdaten mit der Paired-End-Option (z. B. -p). Stellen Sie die heruntergeladene GRCh38 GTF-Annotationsdatei bereit, um die korrekten genomischen Merkmalsgrenzen zu definieren.
4. Berechnen Sie die Fragmente pro Kilobase des Transkripts pro Million abgebildeter Reads (FPKM) für jedes Gen. Verwenden Sie FPKM-Werte ausschließlich für beschreibende Visualisierung, PCA, Heatmap-Darstellung und explorative Ausdruckszusammenfassungen; verwenden Sie keine FPKM-Werte als Eingabematrix für DESeq2-Differentialausdruckstests.

6. Analyse der differentiellen Genexpression

1. Führen Sie eine differentielle Expressionsanalyse zwischen den neu diagnostizierten und rückfälligen Gruppen mit R v3.5.0 und dem DESeq2 R-Paket v1.20.0 durch. Importiere die in Schritt 5.2 generierte Roh-Lesezählmatrix in die R-Umgebung und behalte FPKM-Werte nur für Visualisierung und explorative Analysen.
2. Erstelle das spezialisierte Datensatzobjekt, das vom Analysepaket benötigt wird. Führe den spezifischen Befehl aus (z. B. DESeqDataSetFromMatrix()), um die Zähldatenmatrix mit der entsprechenden Beispielmetadatentabelle zu binden.
3. Definieren Sie die experimentelle Designformel innerhalb des Softwareobjekts. Spezifizieren Sie den klinischen Zustand (neu diagnostiziert versus rückfällig) als primäre Vergleichsvariable (z. B. Design = ~ Zustand). Wenn eine technische Batch-Variable in Schritt 4.3 identifiziert wird und nicht vollständig mit dem klinischen Zustand übereinstimmt, wird sie in die Designformel aufgenommen (z. B. Design = ~ Batch + Bedingung).
4. Führen Sie die Kernfunktion der differentiellen Ausdrucksanalyse aus (z. B. DESeq()). Ermöglichen Sie der Software, automatisch Größenfaktorschätzung, Dispersionsschätzung und negativen binomialen Wald-Test²² durchzuführen.
5. Extrahiere die Ergebnistabelle mithilfe der Ergebnisextraktionsfunktion (z. B. Ergebnisse()). Spezifizieren Sie das Kontrastargument, um den genauen Vergleich zu definieren (rückfällig versus neu diagnostiziert).
6. Passen Sie die resultierenden P-Werte an, um die falsche Entdeckungsrate zu kontrollieren. Nutzen Sie das integrierte Benjamini-und-Hochberg-Verfahren, das automatisch vom Softwarepaket²³ angewendet wird.
7. Filtern Sie die extrahierte Ergebnistabelle, um die signifikanten differenziell exprimierten Gene (DEGs) zu isolieren. Jedes Gen mit einem angepassten P-Wert < 0,05 und einer absoluten logarimatischen_2-fachen Änderung > 1 als signifikant differenziell exprimiert zuzuweisen.

7. Genontologie (GO) Anreicherungsanalyse

1. Führen Sie eine Genontologie-(GO)-Anreicherungsanalyse der identifizierten DEGs mit clusterProfiler v3.8.1 und org. durch. Hs.eg.db v3.6.0. Geben Sie die Liste der Entrez-Gen-IDs ein, die den signifikanten DEGs entsprechen, die in Schritt 6.7 identifiziert wurden.
2. Führen Sie die GO-Anreicherungsfunktion aus (z. B. enrichGO()). Spezifizieren Sie die erforderlichen Parameter, einschließlich der entsprechenden Hintergrund-Organismusdatenbank (z. B. OrgDb = org. Hs.eg.db), der spezifischen Ontologiedomäne (Biologischer Prozess, zelluläre Komponente oder Molekulare Funktion) und der angepasste P-Wert-Grenzwert (0,05).
3. Stellen Sie sicher, dass der Algorithmus während der Anreicherungsberechnung die notwendigen Korrekturen vornimmt. Bestätigen Sie, dass die Software intern die Genlängenverzerrung korrigiert und die P-Werte mit der Benjamini-und-Hochberg-Methode²⁴ anpasst.
4. Betrachten Sie GO-Terme mit einem korrigierten P-Wert kleiner als 0,05 als deutlich angereichert. Erstellen Sie ein Punktdiagramm oder ein Balkendiagramm mit den integrierten Visualisierungsfunktionen des Pakets, um die wichtigsten angereicherten GO-Begriffe darzustellen.

8. Kyoto Encyclopedia of Genes and Genomes (KEGG) Signalweganreicherungsanalyse

1. Nutzen Sie eine umfassende Datenbankressource, die sich dem Verständnis der Funktionen von hochrangigen biologischen Systemen widmet, um dysregulierte Signalwege zu identifizieren. Erstellen Sie die gleiche Liste der wichtigsten DEG Entrez-IDs wie in Schritt 7.1.
2. Führen Sie die KEGG-Anreicherungsfunktion (z. B. enrichKEGG()) innerhalb des funktionalen Annotationspakets R aus.
3. Definieren Sie die kritischen Parameter innerhalb des Funktionsaufrufs. Setze den Organismus-Code streng auf Mensch (z. B. Organismus = 'hat') und definiere die P-Wert-Anpassungsmethode (z. B. pAdjustMethod = 'BH').
4. Extrahiere die statistisch signifikanten KEGG-Signalwege. Filtere die Ausgabe so, dass nur jene Wege erhalten bleiben, die einen korrigierten P-Wert kleiner als 0,05 zeigen.
5. Visualisieren Sie die am stärksten angereicherten KEGG-Pfade. Nutzen Sie die integrierten Plottfunktionen (z. B. Dotplot()), um die statistische Signifikanz und Genzählungen zu jedem Pfad abzubilden.

9. Gen-Satz-Anreicherungsanalyse (GSEA)

1. Bereiten Sie die vorab bewertete Genliste vor, die für die Analyse erforderlich ist. Berechnen Sie die Rangordnungsmetrik für alle exprimierten Gene mit dem vorzeichenmäßigen -log10(P-Wert) multipliziert mit dem Vorzeichen der_log-2-fachen Änderung, die aus der differentiellen Expressionsanalyse abgeleitet ist.
2. Starten Sie eine lokale Installation der Broad Institute GSEA-Software v4.2.3. Geben Sie die neu generierte vorgerankte Genliste in die Software-Schnittstelle^ein.
3. Laden Sie die erforderlichen vordefinierten Gensätze herunter. Erwerben Sie die Gene Ontology (GO) und die Kyoto Encyclopedia of Genes and Genomes (KEGG) Datensätze aus der Molecular Signatures Database (MSigDB, Version 7.5.1)²⁶.
4. Konfigurieren Sie die Softwareparameter, um den statistischen Anreicherungstest durchzuführen. Stellen Sie die Anzahl der Permutationen auf 1.000 ein und wählen Sie den Permutationstyp als 'gene_set' aus.
5. Führen Sie den Analysealgorithmus aus, um festzustellen, ob die vordefinierten Gensätze einen statistisch signifikanten, übereinstimmenden Unterschied zwischen den neu diagnostizierten und rückfälligen biologischen Zuständen zeigen.
6. Bewerten Sie die statistische Signifikans der generierten Anreicherungsprofile. Definiere signifikante Gensätze mit strengen Schwellenwerten: einen normalisierten Anreicherungswert (NES) Absolutwert > 1,0, einen nominalen P-Wert < 0,05 und einen False Discovery Rate (FDR) q-Wert < 0,25.

10. Analyse des Protein-Protein-Interaktionsnetzwerks (PPI)

1. Zugriff auf die STRING-Datenbank für bekannte und vorhergesagte Protein-Protein-Interaktionen. In diesem Workflow wurde die PPI-Analyse mit STRING v11.5²⁷ durchgeführt.
2. Geben Sie die Liste der Entrez-Gen-IDs oder offiziellen Gensymbole für die signifikant differenziell exprimierten Gene (identifiziert in Schritt 6.7) in die Datenbanksuchoberfläche ein. Wähle Homo sapiens als Zielorganismus aus.
3. Konfigurieren Sie die Netzwerkkonstruktionsparameter, um sicherzustellen, dass hochwertige Interaktionen abgerufen werden. Setzen Sie den erforderlichen Mindest-Interaktionsscore auf eine hohe Konfidenzschwelle (Wert > 0,700).
4. Exportiere die resultierenden Netzwerkdaten der Interaktion in ein lokales Verzeichnis. Speichern Sie die Interaktionskarte als Standardtabellendatei (z. B. TSV-Format).
5. Importieren Sie die exportierten Interaktionsdaten in Cytoscape v3.9.1 zur Netzwerkvisualisierung und -analyse²⁸.
6. Filtern Sie das aufgebaute Netzwerk, um die Visualisierungsklarheit zu verbessern, und heben Sie wichtige regulatorische Zentren hervor. Entfernen Sie alle losgelösten Knoten oder Orphan-Gene, die keine kontinuierlichen Interaktionen zeigen, die die festgelegte Konfidenzschwelle erfüllen.

Ergebnisse

Klinische Kohorten- und Sequenzierungsvalidierung

Die erfolgreiche Ausführung des vorgelagerten RNA-Extraktions- und Bibliotheksvorbereitungsprotokolls (Abbildung 1) wurde durch Sequenzierungs-, Ertrags- und Qualitätsmetriken bestätigt. In diesem repräsentativen Datensatz ergaben Knochenmarkproben von fünf neu diagnostizierten AML-Patienten und vier rückfälligen AML-Patienten durchschnittlich etwa 6,0 GB Rohdaten pro Probe. Die Qualitätskontrollbewertung (Tabelle 2) bestätigte, dass die Basisqualität und die Lesetiefe die erforderlichen Schwellenwerte für die nachgelagerte bioinformatische Analyse erfüllten. Niedrige RNA-Integrität (zum Beispiel RIN < 6,0), niedrige Mapping-Raten oder hoher Transkriptdegradations-Bias würden eine suboptimale Eingabequalität darstellen und die Zuverlässigkeit der nachgelagerten differentiellen Expressionsanalyse beeinträchtigen.

Globale transkriptomische Varianz und PCA

Um die globale Transkriptomvarianz zu bewerten und die klinische Gruppierung zu untersuchen, wurde PCA an den normalisierten Expressionsdaten durchgeführt. In diesem repräsentativen Datensatz zeigten die neu diagnostizierten und rückfälligen Gruppen eine Trennung im zweidimensionalen Raum (Abbildung 2A)²⁰, wobei PC1 und PC2 jeweils 23,82 % bzw. 18,75 % der Gesamtvarianz ausmachten. Die Venn-Diagramme in Abbildung 2B,C bieten eine zusätzliche beschreibende Zusammenfassung der Gene, die in Proben der neu diagnostizierten und rückfälligen Gruppen nachgewiesen wurden, und unterstützen Proproduzierbarkeitsprüfungen auf Probenebene vor der nachgelagerten differentiellen Expressionsanalyse. Da die Kohorte klein und ungepaart war, wurde die PCA-Trennung als illustrative Workflow-Ausgabe und nicht als endgültige Belege für krankheitszustandsspezifische Biologie interpretiert.

Analyse des differentiellen Expressionsgens (DEG)

Die Anwendung der festgelegten Protokollschwellenwerte (|log₂FC| ≥ 1 und angepasster P-Wert ≤ 0,05) auf den DESeq2-Ausgang ergab 2.025 DEGs, die 772 hochregulierte und 1.253 herunterregulierte Gene in der rückfälligen Gruppe umfassen (Abbildung 3A). Kandidatentranskripte mit hoher Variation umfassten FOXC1 (log₂FC = 7,55, P = 4,92 x ^10-5), HOXA11 (log₂FC = 7,76), HOXA11-AS (log₂FC = 7,23) und AXL (log₂FC = 3,50) sowie herunterregulierte RHOB, PTX3 und CXCL8. Die vorhandene Literatur verbindet mehrere dieser Gene mit AML-Stammismus, Signalisierung oder Therapie-Response^13,29; Der aktuelle Workflow identifiziert sie jedoch nur als rückfällig-assoziierte Kandidaten-Transkripte. Jede eindeutige mechanistische Rolle bei klinischer Resistenz erfordert eine anschließende unabhängige funktionelle Validierung.

Funktionale und Pathway-Anreicherung (GO, KEGG und GSEA)

Das funktionale Annotationsprotokoll ordnete die DEGs breiteren biologischen Systemen zu. Die GO-Analyse identifizierte die Anreicherung von Termen im Zusammenhang mit der kleinen GTPase-vermittelten Signaltransduktion, dem Metallionentransport und der Chromatinmontage (Abbildung 4A–C). Die KEGG-Wegkartierung identifizierte Zusammenhänge mit ECM-Rezeptor-Interaktionen und Zytokin-Zytokin-Rezeptor-Interaktionen (Abbildung 4D). GSEA zeigte eine Anreicherung von RNA-Biosynthetikprozessen in der rückfälligen Gruppe und die Anreicherung der Energiestoffwechselwege in der neu diagnostizierten Gruppe (Abbildung 5A). Diese Anreicherungsergebnisse bieten eine beschreibende Roadmap veränderter Gensätze und sollten als hypothesenbildende Assoziationen und nicht als nachgewiesene Treiber von Rückfällen interpretiert werden.

Aufbau des Protein-Protein-Interaktionsnetzwerks (PPI)

Das anfängliche STRING-Netzwerk enthielt 56 Knoten und 193 Interaktionen. Nach der Entfernung von losgelösten oder verwaisten Knoten enthielt das dargestellte Cytoscape-Subnetzwerk 42 Knoten und 136 Interaktionen (Abbildung 5B). Die Netzwerkmodulanalyse priorisierte TP53, CCL2, CXCL8 und IL6 als zentrale mathematische Hubs mit der höchsten Anzahl an Interaktionen. Da das PPI-Netzwerk auf datenbankvorhergesagten Interaktionsscores (z. B. ATF3-Score: 0,982) basiert, sollte die Hub-Identifikation als Zielpriorisierung für zukünftige empirische Studien interpretiert werden und nicht als direkte Belege für p53-vermittelte Apoptose-Evasion oder andere Resistenzmechanismen.

Die in diesem Protokoll generierten Rohdaten zur RNA-Sequenzierung wurden im Figshare-Repository abgelegt und sind öffentlich über das folgende DOI zugänglich: https://doi.org/10.6084/m9.figshare.30655814. Die verarbeiteten Daten und die zugehörigen Analysedateien sind im Artikel und den ergänzenden Materialien enthalten. Repräsentative Kommandozeilenparameter und Analyseeinstellungen, die zur Reproduktion des Berechnungsworkflows verwendet werden, werden als Ergänzungsdatei 1 bereitgestellt. Alle Daten, die die Ergebnisse dieser Studie stützen, sind uneingeschränkt verfügbar.

Patientenidentifikation	Alter (Jahre)	Geschlecht	Molekulare Mutationen	Überleben/Nachbeobachtung (Monate)	Klinischer Status
R_AML_1	70	Männlich	FLT3-ITD (+)	22	Verstorben
R_AML_2	29	Weiblich	NPM1 (+)	11	Am Leben
R_AML_3	40	Männlich	CEBPA (+)	17	Am Leben
R_AML_4	55	Weiblich	Dreifach negativ*	24	Verstorben

Tabelle 1: Demografische und klinische Merkmale der Patienten in der rückfälligen AML-Gruppe (R_AML). Tabelle 1 fasst die demografischen und klinischen Merkmale der in der repräsentativen Analyse verwendeten rezidiven AML-Kohorte zusammen, einschließlich klinischer Merkmale auf Patientenebene, die für die Interpretation des transkriptomischen Workflows relevant sind.

Beispiel	Bibliothek	Raw_reads	Raw_bases	Clean_reads	Clean_bases	Error_rate	F20	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

Tabelle 2: Zusammenfassung der Datenqualität. Tabelle 2 zeigt Sequenzierungsqualitätsmetriken für jede Probe, einschließlich Leserbeute, Basisqualität, GC-Gehalt und mappingbezogenen Qualitätskontrollinformationen, die zur Bestimmung der Eignung der Proben für eine nachgelagerte Analyse verwendet werden.

Abbildung 1: Arbeitsablauf des Protokolls. Der Arbeitsablauf fasst die wichtigsten experimentellen und rechnerischen Phasen zusammen, darunter klinische Probenentnahme, RNA-Qualitätskontrolle, Bibliotheksvorbereitung und -sequenzierung, Leseverarbeitung und -ausrichtung, Transkriptquantifizierung, Differentialexpressionsanalyse, GO/KEGG-Anreicherung, GSEA und PPI-Netzwerkaufbau. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzusehen.

Abbildung 2: Quantitative Analyse der Stichproben. (A) Die Hauptkomponentenanalyse (PCA) wurde durchgeführt, um Intergruppenunterschiede und die Reproduzierbarkeit innerhalb der Gruppe zu bewerten. PCA wurde mit linearen algebraischen Methoden durchgeführt, die auf normalisierten Genexpressionswerten über alle Proben basieren. (B, C) Venn-Diagramme, die Gene zeigen, die über Proben in der AML- bzw. R_AML-Gruppe nachgewiesen wurden. Probenbeschränkte Regionen zeigen Gene an, die in einzelnen Proben nachgewiesen wurden, während überlappende Bereiche Gene darstellen, die häufig über zwei oder mehr Proben nachgewiesen werden. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzusehen.

Abbildung 3: Analyse der differenziellen Genexpression. (A) Balkendiagramm, das die Anzahl der differenziell exprimierten Gene (DEGs) zwischen Vergleichsgruppen zeigt, identifiziert durch DESeq2 mit Schwellenwerten des angepassten P-Werts ≤ 0,05 und |log₂FoldChange| ≥ 1. (B) Vulkan-Plot von DEGs. Die x-Achse steht fürlog-2-FoldChange-Werte, und die y-Achse für -log10(P-Wert). Blaue gestrichelte Linien zeigen die Schwellenlinien an, die für die DEG-Auswahl verwendet werden. (C) Hierarchische Clustering-Heatmap der DEGs. Die x-Achse bezeichnet die Stichprobennamen, und die y-Achse zeigt normalisierte Ausdruckswerte der DEGs. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzusehen.

Abbildung 4: Funktionelle Anreicherungsanalyse unterschiedlich exprimierter Gene. (A) GO-Anreicherungsbalken-Plot. Die x-Achse stellt GO-Terme dar, und die y-Achse zeigt die Anreicherungssignifikans, ausgedrückt als -log10(padj). Die Farben stehen für BP (Biologischer Prozess), CC (Zelluläre Komponente) und MF (Molekulare Funktion). (B) Setze die Anreicherungsblasen-Plot. Die x-Achse stellt das Verhältnis der zu jedem GO-Term annotierten DEGs relativ zur Gesamtzahl der DEGs dar, und die y-Achse zeigt die GO-Terme an. Die Blasengröße entspricht der Anzahl der annotierten Gene, und Farbverläufe stellen die Bedeutung der Anreicherung dar. (C) KEGG-Anreicherungsbar-Plot. Die x-Achse stellt KEGG-Signalwege dar, und die y-Achse bezeichnet die Anreicherungsbedeutung. (D) KEGG-Anreicherungsblasen-Plot. Die Blasengröße zeigt die Anzahl der annotierten Gene an, und Farbverläufe spiegeln die Bedeutung der Anreicherung wider. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzusehen.

Abbildung 5: GSEA-Anreicherung und Analyse des Protein-Protein-Interaktionsnetzwerks (PPI). (A) Bar-Plot, das normalisierte Anreicherungswerte (NES) für ausgewählte signifikante Gensätze zeigt. Positive NES-Werte deuten auf Anreicherung in der R_AML-Gruppe hin, während negative NES-Werte auf Anreicherung in der neu diagnostizierten AML-Gruppe hinweisen. (B) Protein-Protein-Interaktionsnetzwerk (PPI). Jeder Knoten stellt ein Protein dar, und jede Kante kennzeichnet eine Wechselwirkung zwischen verbundenen Proteinen. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzusehen.

Diskussion

Kritische Schritte im Protokoll

Die erfolgreiche Ausführung dieses bioinformatischen Workflows beruht auf mehreren kritischen Schritten. Erstens sind das sofortige Einfrieren und die ordnungsgemäße Lyse des Knochenmarkaspirats (Schritt 1.6) von größter Bedeutung, da das Knochenmarkmikroumfeld reich an Ribonukleasen ist, die die transkriptomische Integrität³⁰ schnell abbauen können. Während der Rechenphase ist die korrekte Konfiguration der experimentellen Designformel im DESeq2-Paket (Schritt 6.3) entscheidend für eine genaue differentielle Expression, insbesondere wenn der klinische Zustand (neu diagnostiziert versus rückfällig) gegenübergestellt wird, während potenzielle Störfaktoren berücksichtigt werden. Schließlich ist die Anwendung rigoroser False Discovery Rate (FDR)-Schwellenwerte während der Gene Set Enrichment Analysis (GSEA) (Schritt 9.6) ein entscheidender statistischer Kontrollpunkt, um eine Überinterpretation falsch-positiver funktionaler Netzwerke zu verhindern.

Modifikationen und Fehlerbehebung

Eine häufige Herausforderung bei dieser Methode ist das Vorhandensein von Batch-Effekten, die häufig auftreten, wenn klinische Proben über längere Zeiträume entnommen und sequenziert werden. Chargenvariablen sollten vor der Analyse definiert werden, einschließlich Probenentnahmedatum, RNA-Extraktionsdatum, Bibliotheksvorbereitungscharge, Sequenzierungslane und Sequenzierungslauf. Wenn PCA- oder Stichprobenkorrelationsanalysen eine Clusterbildung basierend auf dem Sequenzierungsdatum oder einer anderen technischen Variable statt klinischen Phänotyps aufzeigen, sollten Nutzer das Protokoll modifizieren, indem sie die Batch-Variable in die Differentialausdrucksdesignformel aufnehmen, wenn es statistisch möglich ist, oder indem sie Batch-Korrekturalgorithmen wie ComBat oder SVA vor Visualisierung³¹ anwenden. Wenn dieses Protokoll auf Vollblutaspiraten anstelle von Knochenmarkaspiraten angewendet wird, ist eine wesentliche Modifikation die Einbeziehung eines Globin-mRNA-Depletionsschritts während der Bibliotheksvorbereitung, um zu verhindern, dass hochhäufige Globintranskripte die Sequenzierungslesetiefe monopolisieren. Softwareversionen und Hauptparameter für den repräsentativen Workflow wurden wie folgt ergänzt: 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, gepaarte 150-bp-Sequenzierung, GSEA v4.2.3 mit 1.000 Gensatz-Permutationen, MSigDB v7.5.1, STRING v11.5 mit hochkonfidenziellen Interaktionen und Cytoscape v3.9.1. Repräsentative Kommandozeilenparameter und Analyseeinstellungen werden in Supplementary File 1 bereitgestellt.

Einschränkungen der Methode

Obwohl umfassend, weist dieses Protokoll inhärente methodische Einschränkungen auf. Erstens verwendet sie Bulk-RNA-Sequenzierung, die das durchschnittliche transkriptomische Profil des gesamten Knochenmarkaspirats erfasst und keine räumliche Einzelzellauflösung bietet. Daher kann der Workflow nicht feststellen, ob eine hochregulierte, rückfallassoziierte Signatur von Leukämiestammzellen, Stromalzellen, Immunzellen oder Veränderungen der Zelltypzusammensetzung³² stammt. Zweitens ist der repräsentative Datensatz klein (n = 9) und ungepaart, was die statistische Robustheit einschränkt und eine eindeutige kausale Schlussfolgerung verhindert. Drittens ist der Workflow rein in Silico. Es erzeugt Kandidat-Regulationszentren und Signalwege, kann deren funktionelle Notwendigkeit bei Chemoresistenz jedoch nicht unabhängig validieren, ohne orthogonale in vitro oder in vivo experimentelle Validierung.

Jüngste Einzelzell- und Einzelzell-genomische Studien haben das AML-Referenzrahmen erweitert, indem sie Zellzustandsheterogenität, klonale Architektur und therapieassoziierte Evolution mit höherer Auflösung^{aufgeklärt haben.} Diese Ansätze ergänzen den hier beschriebenen Bulk-RNA-Seq-Workflow: Bulk-Sequenzierung bietet eine praktische und kosteneffiziente Screening-Strategie für Kohorten-Transkriptom-Signaturen, während Einzelzell- und Multi-Omic-Methoden in Folgestudien verwendet werden können, um Kandidatensignale bestimmten malignen oder mikro-umweltlichen Zellpopulationen zuzuweisen.

Bedeutung im Hinblick auf bestehende Methoden

Trotz dieser Einschränkungen bietet diese transkriptomische Pipeline Vorteile gegenüber alternativen diagnostischen und analytischen Techniken. Traditionelle klinische Beurteilungen des AML-Rückfalls basieren oft auf gezielten multiplexen qPCR-Panels oder Standard-Durchflusszytometrie. Obwohl sie für schnelle Diagnostik nützlich sind, sind diese gezielten Methoden durch vordefinierte Sonden eingeschränkt und können nur bekannte Widerstandsmarker¹⁹ bewerten. Durch die Verwendung unvoreingenommener, genomweiter Transkriptomsequenzierung in Kombination mit Netzwerkanalyse kann dieses Protokoll neuartige Transkripte und systemweite Assoziationen nominieren, die bestehende gezielte Methoden möglicherweise übersehen.

Bedeutung und potenzielle Anwendungen

Die in diesem Protokoll beschriebene Methodik ist relevant für die translationale Hämatologie und personalisierte Medizin, da sie transkriptomische Signaturen, die mit Rückfällen assoziiert sind, für weitere Studien priorisieren kann. Eine mögliche nachgelagerte Anwendung ist die Nomination von Oberflächenantigenen oder Immunumgehungsweges, die während des Rückfalls entstehen. Solche Kandidaten könnten die Gestaltung zukünftiger Validierungsstudien beeinflussen und, falls sie experimentell bestätigt werden, zur Entwicklung von Immuntherapien der nächsten Generation beitragen, einschließlich der CAR-T- oder CAR-NK-Zellstrategien³⁷. Diese translationalen Anwendungen sollten eher als hypothesengenerierend als als etablierte Schlussfolgerungen aus dem vorliegenden Datensatz betrachtet werden.

Materialien

Liste der in diesem Artikel verwendeten Materialien

Name	Unternehmen	Katalognummer	Kommentare
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