Profilazione trascrivomica e analisi bioinformatica di campioni di midollo osseo per identificare le firme di resistenza alla chemioterapia nella leucemia mieloide acuta

Sommario

Questo protocollo presenta un flusso di lavoro bioinformatico standardizzato per analizzare le alterazioni trascridemiche nella leucemia mieloide acuta (LMA). L'obiettivo è confrontare i campioni di midollo osseo appena diagnosticati e quelli recidivati e dare priorità alle firme molecolari associate alla resistenza alla chemioterapia e alla progressione della malattia per ulteriori indagini.

Abstract

La leucemia mieloide acuta (LMA) è una malignità ematologica altamente eterogenea in cui la ricaduta e la chemioresistenza acquisita rimangono cause principali di fallimento del trattamento. Questo articolo presenta un protocollo bioinformatico per l'analisi trascritomica degli aspiratori di midollo osseo. L'obiettivo principale del protocollo è fornire un flusso di lavoro standardizzato per identificare le firme molecolari associate alla progressione della malattia e alla resistenza terapeutica nella LMA recidivante. La pipeline dettaglia le procedure computazionali per confrontare campioni di midollo osseo non accoppiati, dimostrate utilizzando dati di sequenziamento di cinque casi di nuova diagnosi e quattro casi ricaduti. Questo metodo delinea i passaggi essenziali per elaborare i dati di sequenziamento dell'RNA, eseguire l'analisi differenziale dell'espressione genica e condurre valutazioni funzionali a valle. Applicando questo flusso di lavoro, sono stati identificati 2.025 geni differenzialmente espressi (DEG), inclusi FOXC1, HOXA11, HOXA11-AS e AXL, come trascrizioni candidate associate alla ricaduta in questo dataset rappresentativo. Le analisi funzionali e di rete hanno ulteriormente dato priorità ai set genici e ai hub di interazione relativi alla segnalazione di piccole GTPasi, alla segnalazione infiammatoria, alle interazioni della matrice extracellulare e ai processi biosintetici dell'RNA. Nel complesso, questa metodologia fornisce una pipeline computazionale riproducibile per mappare le firme trascrittomiche associate all'AML ricaduta e per generare ipotesi che richiedono una successiva validazione sperimentale.

Introduzione

La leucemia mieloide acuta (LMA) è un gruppo di neoplasie maligne clonali originate da cellule staminali/progenitrici ematopoietiche, caratterizzate da una proliferazione anomala di cellule mieloidi immature nel midollo osseo e dalla soppressione della differenziazione ematopoietica ^1,2,3. Sebbene la chemioterapia di induzione standard attuale (come la citarabina combinata con antracicline) possa indurre una remissione completa nella maggior parte dei pazienti, il tasso di ricaduta rimane alto fino al 50%–70%, e la prognosi per i pazienti recidivanti è significativamente^{compromessa 4,5,6}. La resistenza acquisita alla chemioterapia è associata al fallimento del trattamento AML, rendendo necessaria un'analisi approfondita delle firme molecolari collegate a questo processo per migliorare le strategie terapeutiche e i tassi di sopravvivenza dei pazienti ^7,8.

Nel più ampio corpus della letteratura, studi esistenti indicano che la resistenza alla chemioterapia nella LMA non si limita alla regolazione all'alza delle pompe di efflusso dei farmaci o al metabolismo anomalo dei farmaci, ma è anche collegata al mantenimento della sopravvivenza delle cellule staminali leucemiache (LSC), alla formazione di fenotipi simili alla transizione epiteliale-mesenchimale (EMT) all'interno del settore ematologico e alle interazioni con il microambiente del midollo^{osseo 9. 10,11}. Ad esempio, le popolazioni di LSC mostrano una grande capacità di auto-rinnovamento e quiescenza, associata a una resistenza intrinseca agli agenti chemioterapeutici specifici del ciclo^{cellulare 12}. Inoltre, la regolazione all'alza delle tirosina chinasi recettori, come AXL, è stata associata alla resistenza nell'AML FLT3-ITD+, che si verifica parallelamente all'attivazione delle vie PI3K/AKT e MAPK e a capacità anti-apoptotiche^{potenziate 13,14}. La riprogrammazione metabolica e la rimodellazione epigenetica sono state riconosciute come importanti assi regolatori nella formazione della resistenza. Le evidenze suggeriscono che le cellule AML durante la ricaduta possano adattarsi allo stress ossidativo indotto dalla chemioterapia e al danno al DNA tramite un'attività di fosforilazione ossidativa aumentata (OXPHOS), rapporti NAD⁺/NADH modulati e stati alterati di modifica degli^istoni 15,16,17. Fattori infiammatori nel microambiente del midollo osseo, come IL-6 e CXCL8, sono anch'essi associati alla sopravvivenza degli LSC e alla resistenza alla chemioterapia, spesso in coordinamento con l'attivazione delle vie di segnalazione STAT3/^{NF-κB 11,18}.

Nonostante questi meccanismi riconosciuti, i cambiamenti trascrittomici associati alla ricaduta dell'AML e alla chemioresistenza rimangono incompletamente caratterizzati, soprattutto quando i segnali metabolici, epigenetici e microambientali del midollo osseo vengono valutati insieme nei campioni clinici. Questo flusso di lavoro risponde a questa esigenza dando priorità ai candidati DEG, ai percorsi e alle reti regolatorie associate alla transizione dalla diagnosi iniziale alla ricaduta clinica. L'approccio integra analisi di espressione differenziale, analisi di arricchimento del set genico (GSEA) e costruzione di reti di interazione proteina-proteina (PPI) per mappare la riprogrammazione trascrizionale a livello di sistema e per nominare candidati per una successiva validazione meccanicistica.

L'obiettivo generale di questo metodo è presentare una pipeline bioinformatica standardizzata e riproducibile per confrontare i trascritti in massa di campioni di midollo osseo AML appena diagnosticati e quelli recidivati. La ragione dietro l'uso di questa tecnica in silico è la sua capacità di catturare eventi trascrizionali imparziali e a livello genomico, andando oltre i limiti delle analisi a singolo percorso per dare priorità a complesse reti regolatorie multidimensionali. Questa tecnica offre vantaggi significativi rispetto a metodi alternativi, come i microarray o i pannelli multiplex qPCR mirati, fornendo un intervallo dinamico più elevato, la capacità di rilevare nuovi trascritti e una quantificazione precisa dell'espressione genica senza le limitazioni delle sonde^{predefinite 19,20}. Per determinare se questo metodo sia appropriato per la loro applicazione, i lettori dovrebbero notare che questa pipeline è specificamente progettata per i ricercatori che elaborano dati di sequenziamento di RNA in massa provenienti da coorti cliniche accoppiate o non accoppiate, come le aspirazioni tissutali. È adatto per identificare firme ampie associate alla resistenza e reti regolatorie candidate, mentre i ricercatori che richiedono risoluzione specifica per tipo di cellula o spaziale dovrebbero utilizzare flussi di lavoro complementari di sequenziamento a singola cella o spaziale. In definitiva, questo protocollo computazionale consente la priorità dei geni candidati, delle vie e delle reti regolatorie per ulteriori indagini sperimentali.

Protocollo

Tutti i metodi che prevedevano l'uso di campionamento di tessuto umano sono stati eseguiti in conformità con le linee guida istituzionali e la Dichiarazione di Helsinki (revisionata nel 2013). I campioni clinici di midollo osseo sono stati ottenuti con l'approvazione del comitato etico istituzionale (n. di approvazione TY-ZKY2024-116-01 e TY-ZKY2024-116-02).

1. Raccolta clinica di campioni e classificazione dei pazienti

1. Seleziona campioni di aspirazione del midollo osseo da pazienti formalmente diagnosticati con leucemia mieloide acuta (LMA) basandosi sui criteri di classificazione dell'Organizzazione Mondiale della Sanità (OMS).
2. Applicare criteri specifici di inclusione ed esclusione durante la selezione dei pazienti per garantire l'omogeneità e la riproducibilità delle coorti. Includere i pazienti adulti con LMA primaria ed escludere i pazienti con LMA secondaria, leucemia promielocita acuta o una storia precedente di altre neoplasie (Tabella 1).
3. Assegnare il campione raccolto al gruppo di nuova diagnosi se il paziente presenta AML non trattato al momento della diagnosi clinica iniziale.
4. Assegnare il campione raccolto al gruppo recidivo se il paziente mostra una ricomparsa di blasti leucemici nel sangue periferico o blasti superiori al 5% nel midollo osseo dopo una remissione completa documentata.
5. Raccogliere i campioni residui di midollo osseo deidentificati immediatamente dopo la procedura clinica di aspirazione di routine del midollo osseo.
   NOTA: Nella dimostrazione di questo protocollo specifico, sono stati raccolti nove campioni consecutivi (cinque appena diagnosticati e quattro ricaduti) da settembre 2024 a settembre 2025. Poiché questo utilizzava solo campioni clinici residui de-identificativi, il comitato etico ha revocato il requisito del consenso informato scritto.
6. Processare immediatamente l'aspirazione del midollo osseo raccolto per la conservazione dell'RNA utilizzando un metodo standard di lisi fenolico-guanidinio²¹.
   1. Trasferire il midollo osseo aspirato appena raccolto in un tubo di raccolta contenente un anticoagulante. Scuoti il tubo per mescolare accuratamente aspirato e anticoagulante.
   2. Estrarre un volume misurato del campione anticoagulato e aggiungerlo a un reagente commerciale per la lisi fenolico-guanidinio. Mantenere un rapporto di volume di 3 parti di reagente di lisi per 1 parte di campione.
      ATTENZIONE: Il reagente per la lisi fenol-guanidinio contiene sostanze chimiche tossiche e corrosive che possono causare gravi ustioni e danni ai tessuti. Esegui tutta la gestione dei reagenti all'interno di una cappa chimica indossando l'attrezzatura di protezione individuale appropriata.
   3. Scuoti vigorosamente il tubo per omogeneizzare completamente il campione e il reagente per la lisi. Assicurati che la miscela sia completamente frullata e verifica che non rimangano coaguli visibili nella soluzione.
   4. Congelare immediatamente il campione omogeneizzato immergendo il tubo nell'azoto liquido.
      ATTENZIONE: L'azoto liquido è estremamente freddo e può causare gravi congelamenti al contatto. Indossa guanti criogenici e una visiera completa quando manipoli azoto liquido.
7. Trasferire i campioni congelati rapidamente in un congelatore a -80 °C per la conservazione a lungo termine prima della pipeline di isolamento e sequenziamento del trascrittoma a valle dell'RNA. Questo rappresenta un punto sicuro in cui l'esperimento può essere messo in pausa e riavviato successivamente.
   NOTA: Il flusso di lavoro presentato in questo protocollo si concentra interamente sulla generazione di firme di resistenza computazionale. Non è stata effettuata alcuna validazione sperimentale indipendente, come la PCR quantitativa in tempo reale (RT-qPCR), sui geni chiave espressi differenzialmente identificati attraverso questa specifica pipeline.

2. Controllo qualità dell'RNA e preparazione della biblioteca

1. Valutare l'integrità dell'RNA utilizzando un sistema di elettroforesi capillare microfluidica. Per questo flusso di lavoro rappresentativo, includere campioni di RNA con un numero di integrità dell'RNA (RIN) ≥ 6,0, un rapporto A260/280 tra 1,8 e 2,1, e nessun picco di degradazione visibile. Registra i rapporti misurati di RIN e purezza per ogni campione prima della preparazione della biblioteca.
2. Input di 1 μg di RNA totale per campione per la preparazione della biblioteca. Purificare l'mRNA dall'RNA totale utilizzando sfere magnetiche attaccate a oligo-T per arricchire i trascriti con coda di poliA.
3. Frammenta l'mRNA arricchito usando cationi divalenti. Incubare la miscela a 94 °C per 15 minuti in un buffer di reazione di sintesi a primo filamento 5X.
4. Sintetizza il cDNA del primo filamento usando primer esameri casuali e una trascridasi inversa priva di attività RNasi H.
5. Degradare il filamento modello di RNA usando RNase H. Sintetizzare il cDNA di secondo filamento utilizzando DNA polimerasi I e dNTP in un sistema di reazione da 20 μL.
6. Incubare la reazione di sintesi del secondo filamento a 16 °C per 1 ora. Centrifuga brevemente la miscela di reazione a 2.000 x g per raccogliere il liquido sul fondo del tubo.
7. Convertire le sporgenti rimanenti in estremità smussate utilizzando attività di esonucleasi e polimerasi. Adenilate le estremità di 3' dei frammenti di DNA e legare gli adattatori con strutture a forno per prepararsi all'ibridazione.
8. Purificare i frammenti della libreria utilizzando perle di immobilizzazione magnetiche in fase solida reversibile per selezionare preferenzialmente frammenti di cDNA lunghi 370–420 bp.
9. Esegui lavaggi all'etanolo durante la purificazione delle perline. Centrifugare i tubi a 2.000 x g per 30 secondi per raccogliere e rimuovere completamente qualsiasi etanolo residuo prima dell'eluzione finale.
10. Esegui amplificazione PCR utilizzando una DNA polimerasi ad alta fedeltà, primer PCR universali e primer indice specifici per campione.
11. Eseguire il profilo termico PCR con una denaturazione iniziale a 98 °C per 30 secondi. Segui con 12 cicli a 98 °C per 10 secondi, 60 °C per 30 secondi e 72 °C per 30 secondi, concludendo con un'estensione finale a 72 °C per 5 minuti.
12. Purifica di nuovo i prodotti PCR usando le sfere magnetiche. Applicare gli stessi parametri di centrifugazione del passaggio 2.9 per ottenere la libreria finale.
13. Quantifica la concentrazione iniziale della libreria usando un fluorometro. Diluire la libreria finale a una concentrazione di 1,5 ng/μL.
14. Mescola accuratamente la libreria diluita. Centrifugare la miscela a 10.000 x g per 1 minuto a 4 °C per rimuovere eventuali residui prima dell'analisi finale.
15. Valuta la dimensione dell'inserto della libreria utilizzando il sistema di elettroforesi capillare microfluidica.
16. Quantificare accuratamente la concentrazione effettiva della libreria tramite PCR quantitativa in tempo reale (qRT-PCR) dopo aver confermato che la dimensione dell'inserto soddisfi le aspettative. Assicurarsi che la concentrazione sia superiore a 1,5 nM per garantire stabilità della libreria e qualità di sequenziamento.
    NOTA: Questo rappresenta un punto sicuro in cui l'esperimento può essere messo in pausa. Le librerie preparate possono essere memorizzate a -20 °C fino al clustering e al sequenziamento.

3. Clustering e sequenziamento del trascrittoma

1. Effettuare il clustering dei campioni codificati per indice su un sistema automatizzato di generazione di cluster. Usa un kit di cluster a estremità abbinata commerciale secondo le istruzioni del produttore.
2. Sequenziare le preparazioni della libreria su una piattaforma di sequenziamento ad alto throughput dopo la generazione di cluster con successo. Genera 150 letture a coppie di base (bp).

4. Controllo della qualità dei dati e mappatura delle letture

1. Valuta la qualità dei dati grezzi (formato FASTQ) utilizzando fastp v0.23.2 per il controllo qualità e il filtraggio della lettura grezza. Registra i parametri della riga di comando in un registro di analisi. In questo workflow, le letture pulite venivano generate rimuovendo le letture contenenti adattatore, le letture contenenti sequenze poly-N e le letture di bassa qualità usando impostazioni di filtraggio identiche tra i campioni. Un comando di fine accoppiata rappresentativa è fornito nel File Supplementare 1.
2. Elabora le letture grezze tramite un software di pre-elaborazione automatica. Ottieni letture pulite rimuovendo le letture contenenti adattatori, quelle contenenti sequenze poly-N e le letture di bassa qualità. Utilizzare parametri di filtraggio identici per tutti i campioni e registrare il numero di lettura trattenuto, Q20, Q30 e contenuto GC dopo il filtraggio.
3. Calcola il contenuto di Q20, Q30 e GC dei dati puliti. Definire le potenziali variabili batch prima dell'analisi a valle, inclusi data di raccolta del campione, data di estrazione RNA, batch di preparazione della libreria, corsia di sequenziamento e esecuzione di sequenziamento.
4. Valutare gli effetti dei lotti tramite PCA e analisi di correlazione campione-campione utilizzando valori di espressione normalizzati. Se i campioni raggruppano principalmente per variabili tecniche piuttosto che per stato clinico, documentare la variabile interessata e includerla come covariata nella formula di progettazione dell'espressione differenziale oppure applicare un metodo di aggiustamento batch consolidato prima della visualizzazione a valle.
5. Acquisire il genoma di riferimento (Homo sapiens, GRCh38) e i corrispondenti file di annotazione genica Ensembl release 109 per l'allineamento della lettura.
6. Costruisci l'indice del genoma di riferimento usando HISAT2 v2.0.5.
7. Allinea le letture pulite a estremità accoppiata al genoma di riferimento usando HISAT2 v2.0.5. Utilizzare questo approccio di allineamento consapevole dello splice-aware per generare un database di giunzioni di splice basato sul file di annotazione del modello genico.

5. Previsione di trascrizioni innovative e quantificazione dell'espressione genica

1. Assemblare le letture mappate di ogni campione usando StringTie v1.3.3b in un approccio basato sul riferimento. Utilizzare questo strumento per assemblare e quantificare trascritti a lunghezza intera che rappresentano molteplici varianti di splice per ogni locus genico.
2. Contare il numero di letture mappate a ciascun gene usando featureCounts v1.5.0-p3. Utilizzare la matrice di conteggio degli interi grezzi risultante come input per l'analisi differenziale delle espressioni a valle.
3. Configura featureCounts v1.5.0-p3 per i dati di sequenziamento a estremità accoppiata usando l'opzione di estremità accoppiata (ad esempio, -p). Fornire il file di annotazione GTF GRCh38 scaricato per definire i corretti confini delle caratteristiche genomiche.
4. Calcola i Frammenti per Kilobase di trascrizione per milione di letture mappate (FPKM) per ogni gene. Usa i valori FPKM solo per la visualizzazione descrittiva, PCA, visualizzazione della heatmap e riassunti esplorativi delle espressioni; non utilizzare i valori FPKM come matrice di ingresso per il test di espressione differenziale DESeq2.

6. Analisi differenziale dell'espressione genica

1. Eseguire analisi dell'espressione differenziale tra i gruppi appena diagnosticati e quelli recidivati utilizzando R v3.5.0 e il pacchetto DESeq2 R v1.20.0. Importare la matrice grezza del conteggio delle letture generata nel passo 5.2 nell'ambiente R e conservare i valori FPKM solo per visualizzazione e analisi esplorative.
2. Costruisci l'oggetto di dataset specializzato richiesto dal pacchetto di analisi. Esegui il comando specifico (ad esempio, DESeqDataSetFromMatrix()) per associare la matrice dei dati di conteggio con la tabella dei metadati di esempio corrispondente.
3. Definisci la formula di progettazione sperimentale all'interno dell'oggetto software. Specificare lo stato clinico (appena diagnosticato versus recidiva) come variabile principale per il confronto (ad esempio, disegno = ~ condizione). Se una variabile tecnica batch viene identificata nel passo 4.3 e non è completamente confusa con lo stato clinico, includela nella formula di progettazione (ad esempio, design = ~ batch + condizione).
4. Eseguire la funzione centrale di analisi differenziale di espressione (ad esempio, DESeq()). Consentire al software di eseguire automaticamente la stima del fattore dimensione, la stima della dispersione e il test di Wald binomiale negativo²².
5. Estrai la tabella dei risultati utilizzando la funzione di estrazione dei risultati (ad esempio, risultati()). Specifica l'argomento del contrasto per definire il confronto esatto (ricaduta vs nuova diagnosi).
6. Regolare i valori P risultanti per controllare il tasso di false scoperte. Utilizzare la procedura integrata di Benjamini e Hochberg applicata automaticamente dal pacchetto^{software 23}.
7. Filtrare la tabella dei risultati estratti per isolare i geni differenzialmente espressi (DEG) significativi. Assegnare a qualsiasi gene con un valore P aggiustato < 0,05 e un logaritomimetro_{assoluto di variazione di 2} volte > 1 come espresso in modo significativamente differenziale.

7. Analisi dell'arricchimento dell'ontologia genica (GO)

1. Eseguire un'analisi di arricchimento Gene Ontology (GO) dei DEG identificati utilizzando clusterProfiler v3.8.1 e org. Hs.eg.db v3.6.0. Inserisci l'elenco degli ID geni Entrez corrispondenti ai DEG significativi identificati nel passo 6.7.
2. Esegui la funzione di arricchimento GO (ad esempio, enrichGO()). Specificare i parametri richiesti, incluso il database di background degli organismi appropriato (ad esempio, OrgDb = org. Hs.eg.db), il dominio specifico dell'ontologia (Processo Biologico, Componente Cellulare o Funzione Molecolare) e il cutoff del valore P aggiustato (0,05).
3. Assicurarsi che l'algoritmo applichi le correzioni necessarie durante il calcolo dell'arricchimento. Confermare che il software corregga internamente il bias di lunghezza genica e aggiusta i valori P usando il metodo^{Benjamini e Hochberg 24}.
4. Si considerino i termini GO con un valore P corretto inferiore a 0,05 come significativamente arricchiti. Genera un dot plot o un diagramma a barre utilizzando le funzioni di visualizzazione integrate del pacchetto per visualizzare i termini GO arricchiti in cima.

8. Analisi dell'arricchimento delle vie dell'Enciclopedia di Geni e Genomi di Kyoto (KEGG)

1. Utilizzare una risorsa di database completa dedicata alla comprensione delle funzioni di alto livello del sistema biologico per identificare i percorsi disregolati. Prepara la stessa lista di ID DEG Entrez principali utilizzati nel passo 7.1.
2. Esegui la funzione di arricchimento KEGG (ad esempio, enrichKEGG()) all'interno del pacchetto di annotazione funzionale R.
3. Definire i parametri critici all'interno della chiamata di funzione. Imposta il codice dell'organismo strettamente su umano (ad esempio, organismo = 'ha') e definisci il metodo di aggiustamento dei valori P (ad esempio, pAdajustMethod = 'BH').
4. Estrai i percorsi KEGG statisticamente significativi. Filtra l'output per mantenere solo quei percorsi che dimostrano un valore P corretto inferiore a 0,05.
5. Visualizza i percorsi KEGG più arricchiti. Utilizzare le funzioni di graficazione integrate (ad esempio, dotplot()) per mappare la significanza statistica e il conteggio genico associati a ciascun percorso.

9. Analisi dell'arricchimento degli insiemi genici (GSEA)

1. Prepara la lista genetica pre-classificata richiesta per l'analisi. Calcola la metrica di classificazione per tutti i geni espressi usando il segno -log10(P-valore) moltiplicato per il segno della variazione log₂ volte derivata dall'analisi di espressione differenziale.
2. Avvia un'installazione locale del software GSEA del Broad Institute v4.2.3. Inserire la nuova lista genica pre-classificata nella nuova interfaccia^{software 25}.
3. Scarica i set genici predefiniti richiesti. Acquisire i dataset Gene Ontology (GO) e Kyoto Encyclopedia of Genes and Genomes (KEGG) dal Molecular Signatures Database (MSigDB, versione 7.5.1⁾²⁶.
4. Configura i parametri del software per eseguire il test di arricchimento statistico. Imposta il numero di permutazioni a 1.000 e seleziona il tipo di permutazione come 'gene_set'.
5. Eseguire l'algoritmo di analisi per determinare se gli insiemi genici predefiniti mostrano una differenza statisticamente significativa e concordante tra gli stati biologici appena diagnosticati e quelli ricaduti.
6. Valutare la significanza statistica dei profili di arricchimento generati. Definisci insiemi genici significativi usando soglie rigorose: un valore assoluto di punteggio di arricchimento normalizzato (NES) > 1,0, un valore nominale P < 0,05 e un valore q del tasso di falsa scoperta (FDR) < 0,25.

10. Analisi della rete di interazione proteina-proteina (PPI)

1. Accedi al database STRING per le interazioni proteina-proteina note e previste. In questo workflow, l'analisi PPI è stata eseguita con STRING v11.5²⁷.
2. Inserisci l'elenco degli ID geni Entrez o dei simboli genici ufficiali per i geni espressi differenzialmente significativi (identificati nel passaggio 6.7) nell'interfaccia di ricerca del database. Seleziona Homo sapiens come organismo bersaglio.
3. Configurare i parametri di costruzione della rete per garantire che vengano recuperate interazioni di alta qualità. Imposta il punteggio minimo richiesto per l'interazione a una soglia di confidenza elevata (punteggio > 0,700).
4. Esporta i dati di rete di interazione risultanti in una directory locale. Salva la mappa di interazione come file tabulare standard (ad esempio, formato TSV).
5. Importare i dati di interazione esportati in Cytoscape v3.9.1 per la visualizzazione e l'analisi della rete²⁸.
6. Filtra la rete costruita per migliorare la chiarezza della visualizzazione e evidenziare i principali hub regolatori. Rimuovere qualsiasi nodo disconnesso o geni orfani che non mostrano interazioni continue che soddisfino la soglia di fiducia stabilita.

Risultati

Validazione della Coorte Clinica e del Sequenziamento

L'esecuzione con successo del protocollo di estrazione e preparazione della libreria di RNA a monte (Figura 1) è stata confermata dal sequenziamento delle metriche di resa e qualità. In questo dataset rappresentativo, i campioni di midollo osseo di cinque pazienti recentemente diagnosticati con LMA e quattro pazienti ricaduti con LMA hanno fornito in media circa 6,0 GB di dati grezzi per campione. La valutazione del controllo della qualità (Tabella 2) ha confermato che la qualità di base e la profondità di lettura rispettavano le soglie richieste per l'analisi bioinformatica a valle9. Bassa integrità dell'RNA (ad esempio, RIN < 6,0), bassi tassi di mappatura o un alto bias di degradazione dei trascriti rappresenterebbero una qualità di input subottimale e comprometterebbero l'affidabilità dell'analisi differenziale dell'espressione a valle.

Varianza Trascritomica Globale e PCA

Per valutare la varianza globale del trascrittoma e ispezionare il raggruppamento clinico, è stata eseguita PCA sui dati di espressione normalizzata. In questo dataset rappresentativo, i gruppi di nuova diagnosi e ricadute hanno mostrato separazione nello spazio bidimensionale (Figura 2A⁾²⁰, con PC1 e PC2 che rappresentano rispettivamente il 23,82% e il 18,75% della varianza totale. I diagrammi di Venn nella Figura 2B,C forniscono un ulteriore riassunto descrittivo dei geni rilevati tra campioni all'interno dei gruppi di nuova diagnosi e recidivante, supportando i controlli di riproducibilità a livello di campione prima dell'analisi differenziale a valle. Poiché la coorte era piccola e non accoppiata, la separazione delle PCA è stata interpretata come un output illustrativo del flusso di lavoro piuttosto che come una prova definitiva di una biologia specifica dello stato della malattia.

Analisi del Gene di Espressione Differenziale (DEG)

Applicando le soglie stabilite del protocollo (log₂FC| ≥ 1 e valore P aggiustato ≤ 0,05) all'output DESeq2 sono stati individuati 2.025 DEG, comprendenti 772 geni al rialzo e 1.253 al ribasso nel gruppo ricaduto (Figura 3A). I certificati dei candidati con alta variazione includevano FOXC1 (log₂FC = 7,55, P = 4,92 x ^10-5), HOXA11 (log₂FC = 7,76), HOXA11-AS (log₂FC = 7,23) e AXL (log₂FC = 3,50), insieme a RHOB, PTX3 e CXCL8 sotto regolazione. La letteratura esistente collega diversi di questi geni alla stemminalità AML, alla segnalazione o alla risposta^{terapeutica 13,29}; tuttavia, l'attuale flusso di lavoro li identifica solo come trascrizioni di candidati associate a ricadute. Qualsiasi ruolo meccanicistico definitivo nella resistenza clinica richiede una successiva validazione funzionale indipendente.

Arricchimento funzionale e dei percorsi (GO, KEGG e GSEA)

Il protocollo di annotazione funzionale mappava i DEG a sistemi biologici più ampi. L'analisi GO ha identificato un arricchimento di termini relativi alla trasduzione di segnali mediata da piccole GTPasi, al trasporto di ioni metallici e all'assemblaggio della cromatina (Figura 4A–C). La mappatura delle vie KEGG ha identificato associazioni tra interazioni ECM-recettore e citochine-citochine (Figura 4D). GSEA ha mostrato l'arricchimento dei processi biosintetici dell'RNA nel gruppo recidivo e l'arricchimento delle vie metaboliche energetiche nel gruppo appena diagnosticato (Figura 5A). Questi risultati di arricchimento forniscono una tabella di marcia descrittiva dei set genici alterati e dovrebbero essere interpretati come associazioni che generano ipotesi piuttosto che come fattori comprovati della ricaduta.

Costruzione di reti di interazione proteina-proteina (PPI)

La rete STRING iniziale conteneva 56 nodi e 193 interazioni. Dopo la rimozione dei nodi disconnessi o orfani, la sottorete Cytoscape visualizzata conteneva 42 nodi e 136 interazioni (Figura 5B). L'analisi modulare di rete ha dato priorità a TP53, CCL2, CXCL8 e IL6 come centri matematici centrali con il maggior numero di interazioni. Poiché la rete PPI si basa su punteggi di interazione previsti dal database (ad esempio, punteggio ATF3: 0,982), l'identificazione dell'hub dovrebbe essere interpretata come priorità dei target per studi empirici futuri piuttosto che come evidenza diretta di evasione dell'apoptosi mediata da p53 o altri meccanismi di resistenza.

I dati grezzi di sequenziamento dell'RNA generati in questo protocollo sono stati depositati nel repository Figshare e sono accessibili pubblicamente tramite il seguente DOI: https://doi.org/10.6084/m9.figshare.30655814. I dati elaborati e i relativi file di analisi sono inclusi nell'articolo e nei suoi materiali supplementari. I parametri rappresentativi della riga di comando e le impostazioni di analisi utilizzati per riprodurre il flusso di lavoro computazionale sono forniti come Supplementary File 1. Tutti i dati a supporto dei risultati di questo studio sono disponibili senza restrizioni.

ID paziente	Età (Anni)	Sesso	Mutazioni molecolari	Sopravvivenza/Follow-up (mesi)	Stato clinico
R_AML_1	70	Maschio	FLT3-ITD (+)	22	Deceduta
R_AML_2	29	Femmina	NPM1 (+)	11	Vivo
R_AML_3	40	Maschio	CEBPA (+)	17	Vivo
R_AML_4	55	Femmina	Triplo Negativo*	24	Deceduta

Tabella 1: Caratteristiche demografiche e cliniche dei pazienti nel gruppo AML recidivante (R_AML). La Tabella 1 riassume le caratteristiche demografiche e cliniche della coorte AML ricaduta utilizzata nell'analisi rappresentativa, incluse caratteristiche cliniche a livello di paziente rilevanti per l'interpretazione del flusso di lavoro trascritomico.

Esempio	Biblioteca	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

Tabella 2: Riepilogo della qualità dei dati. La Tabella 2 riporta le metriche di qualità del sequenziamento per ogni campione, inclusi rendimento di lettura, qualità di base, contenuto GC e informazioni di controllo qualità relative alla mappatura utilizzate per determinare se i campioni fossero adatti all'analisi a valle.

Figura 1: Flusso di lavoro del protocollo. Il flusso di lavoro riassume le principali fasi sperimentali e computazionali, inclusi la raccolta dei campioni clinici, il controllo qualità dell'RNA, la preparazione e il sequenziamento delle librerie, l'elaborazione e l'allineamento delle letture, la quantificazione dei trascrivi, l'analisi differenziale dell'espressione, l'arricchimento GO/KEGG, la GSEA e la costruzione della rete PPI. Clicca qui per visualizzare una versione più grande di questa figura.

Figura 2: Analisi quantitativa dei campioni. (A) È stata effettuata un'analisi a componenti principali (PCA) per valutare le differenze tra gruppi e la riproducibilità all'interno del campione. La PCA è stata condotta utilizzando metodi algebrici lineari basati su valori di espressione genica normalizzati su tutti i campioni. (B, C) Diagrammi di Venn che mostrano i geni rilevati tra i campioni rispettivamente nei gruppi AML e R_AML. Le regioni a vincoli di campione indicano geni rilevati in singoli campioni, mentre le aree sovrapposte rappresentano geni comunemente rilevati in due o più campioni. Clicca qui per visualizzare una versione più grande di questa figura.

Figura 3: Analisi differenziale dell'espressione genica. (A) Grafico a barre che mostra il numero di geni differenzialmente espressi (DEG) tra i gruppi di confronto, identificati da DESeq2 con soglie di valore P aggiustato ≤ 0,05 e |_{log 2}FoldChange| ≥ 1. (B) Pianta vulcanica dei DEG. L'asse x rappresenta i valori log₂di FoldChange, e l'asse y rappresenta -log10(valore P). Le linee tratteggiate blu indicano le linee di soglia utilizzate per la selezione del DEG. (C) Mappa di calore gerarchica a clustering dei DEG. L'asse x indica i nomi dei campioni, mentre l'asse y mostra valori di espressione normalizzati dei DEG. Clicca qui per visualizzare una versione più grande di questa figura.

Figura 4: Analisi dell'arricchimento funzionale di geni espressi differenzialmente. (A) Trama del bar di arricchimento GO. L'asse x rappresenta i termini GO, e l'asse y mostra la significanza di arricchimento, espressa come -log10(padj). I colori rappresentano BP (Processo Biologico), CC (Componente Cellulare) e MF (Funzione Molecolare). (B) Trama della bolla di arricchimento GO. L'asse x rappresenta il rapporto dei DEG annotati a ciascun termine GO rispetto al numero totale di DEG, mentre l'asse y indica i termini GO. La dimensione della bolla corrisponde al numero di geni annotati, e i gradienti di colore rappresentano la significanza dell'arricchimento. (C) Trama del bar di arricchimento KEGG. L'asse x rappresenta i percorsi KEGG, mentre l'asse y indica il significato per l'arricchimento. (D) Tramo bolla di arricchimento KEGG. La dimensione della bolla indica il numero di geni annotati, e le sfumate di colore riflettono il significato dell'arricchimento. Clicca qui per visualizzare una versione più grande di questa figura.

Figura 5: Arricchimento GSEA e analisi della rete di interazione proteina–proteina (PPI). (A) Grafico a barre che mostra punteggi di arricchimento normalizzato (NES) per set genici significativi selezionati. I valori positivi del NES indicano arricchimento nel gruppo R_AML, mentre valori negativi del NES indicano l'arricchimento nel gruppo con AML appena diagnosticato. (B) Rete di interazione proteina–proteina (PPI). Ogni nodo rappresenta una proteina, e ogni arco indica un'interazione tra proteine connesse. Clicca qui per visualizzare una versione più grande di questa figura.

Discussione

Fasi critiche del protocollo

L'esecuzione efficace di questo flusso di lavoro bioinformatico si basa su diversi passaggi critici. Innanzitutto, il congelamento immediato e la corretta lisi dell'aspirato del midollo osseo (Passo 1.6) sono fondamentali, poiché il microambiente del midollo osseo è ricco di ribonucleasi che possono rapidamente degradare l'integrità^{trascrittomica 30}. Durante la fase computazionale, la corretta configurazione della formula del disegno sperimentale all'interno del package DESeq2 (Passo 6.3) è fondamentale per un'espressione differenziale accurata, soprattutto quando lo stato clinico (appena diagnosticato vs recidivo) viene messo a confronto mentre si considerano le potenziali variabili di confusione. Infine, l'applicazione di rigorose soglie di tasso di false discovery (FDR) durante l'Analisi di Arricchimento dei Setmi Genici (GSEA) (Passo 9.6) è un punto di controllo statistico critico per prevenire l'iper-interpretazione delle reti funzionali falsamente positive.

Modifiche e risoluzione dei problemi

Una sfida comune in questo metodo è la presenza di effetti batch, che si verificano frequentemente quando i campioni clinici vengono raccolti e sequenziati su lunghi intervalli di tempo. Le variabili batch devono essere definite prima dell'analisi, inclusi data di raccolta del campione, data di estrazione RNA, batch di preparazione della libreria, corsia di sequenziamento e esecuzione di sequenziamento. Se l'analisi PCA o di correlazione campionaria rivela un raggruppamento basato sulla data di sequenziamento o su un'altra variabile tecnica piuttosto che sul fenotipo clinico, gli utenti dovrebbero modificare il protocollo includendo la variabile batch nella formula di progettazione dell'espressione differenziale quando statisticamente fattibile o applicando algoritmi di correzione batch, come ComBat o SVA, prima della visualizzazione³¹. Se questo protocollo viene applicato al sangue intero invece che agli aspiratori del midollo osseo, una modifica essenziale è l'inclusione di una fase di esaurimento dell'mRNA della globina durante la preparazione della libreria per impedire che i trascritti di globina ad alta abbondanza monopolizzino la profondità di lettura del sequenziamento. Le versioni software e i parametri principali per il flusso di lavoro rappresentativo sono stati integrati come segue: 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, sequenziamento a estremità accoppiata a 150 bp, GSEA v4.2.3 con 1.000 permutazioni di set genici, MSigDB v7.5.1, STRING v11.5 con interazioni ad alta fiducia e Cytoscape v3.9.1. I parametri rappresentativi della riga di comando e le impostazioni di analisi sono forniti nel File Supplementare 1.

Limitazioni del metodo

Pur essendo completo, questo protocollo presenta limitazioni metodologiche intrinseche. Innanzitutto, utilizza il sequenziamento di RNA in massa, che cattura il profilo trascritomico medio dell'intero midollo aspirato e manca di risoluzione spaziale a singola cellula. Pertanto, il flusso di lavoro non può determinare se una firma associata alla ricaduta aumentata abbia origine da cellule staminali leucemiche, cellule stromali, cellule immunitarie o cambiamenti nella composizione^{cellulare 32}. In secondo luogo, il dataset rappresentativo è piccolo (n = 9) e non accoppiato, il che limita la robustezza statistica e impedisce un'inferenza causale definitiva. Terzo, il flusso di lavoro è puramente in silico. Genera nodi regolatori e vie di segnalazione candidati, ma non può convalidarne indipendentemente la necessità funzionale nella chemioresistenza senza una validazione sperimentale ortogonale in vitro o in vivo.

Studi genomici recenti su singola e singola cellula hanno ampliato il quadro di riferimento AML risolvendo l'eterogeneità cellulare, l'architettura clonale e l'evoluzione associata alla terapia con risoluzioni^{più elevate 33,34,35,36}. Questi approcci sono complementari al flusso di lavoro bulk RNA-seq descritto qui: il sequenziamento in massa fornisce una strategia di screening pratica ed economica per le firme trascridemiche a livello di coorte, mentre metodi a singola cellula e multi-omiche possono essere utilizzati in studi di follow-up per assegnare segnali candidati a specifiche popolazioni cellulari maligne o microambientali.

Importanza rispetto ai metodi esistenti

Nonostante queste limitazioni, questa pipeline trascrittomica offre vantaggi rispetto a tecniche diagnostiche e analitiche alternative. Le valutazioni cliniche tradizionali della ricaduta AML spesso si basano su pannelli multiplex qPCR mirati o citometria a flusso standard. Sebbene utili per diagnostiche rapide, questi metodi mirati sono vincolati da sonde predefinite e possono valutare solo i marcatori di resistenza^{noti 19}. Utilizzando il sequenziamento trascrimico non imparziale e genomico a livello genomico, abbinato all'analisi di rete, questo protocollo può nominare nuovi trascrizioni e associazioni a livello di sistema che i metodi mirati esistenti potrebbero trascurare.

Importanza e potenziali applicazioni

La metodologia descritta in questo protocollo è rilevante per l'ematologia traslazionale e la medicina personalizzata perché può dare priorità alle firme trasscrittomiche associate alla ricaduta per ulteriori studi. Una possibile applicazione a valle è la nomina di antigeni superficiali o vie di evasione immunitaria che emergono durante la ricaduta. Tali candidati potrebbero informare la progettazione di futuri studi di validazione e, se confermati sperimentalmente, potrebbero contribuire allo sviluppo di immunoterapie di nuova generazione, incluse le strategie per le cellule CAR-T^{o CAR-NK 37}. Queste applicazioni traslazionali dovrebbero essere considerate generanti ipotesi piuttosto che conclusioni consolidate dal dataset attuale.

Materiali

Elenco dei materiali utilizzati in questo articolo

Nome	Azienda	Numero di catalogo	Commenti
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