Profilage transcriptomique et analyse bioinformatique d’échantillons de moelle osseuse pour identifier les signatures de résistance à la chimiothérapie dans la leucémie myéloïde aiguë

Résumé

Ce protocole présente un flux de travail bioinformatique standardisé pour analyser les altérations transcriptomiques dans la leucémie myéloïde aiguë (LMA). L’objectif est de comparer les échantillons de moelle osseuse nouvellement diagnostiqués et rechutés, et de prioriser les signatures moléculaires associées à la résistance à la chimiothérapie et à la progression de la maladie pour des investigations ultérieures.

Résumé

La leucémie myéloïde aiguë (LMA) est une malignité hématologique très hétérogène dans laquelle la rechute et la chimiorésistance acquise restent des causes majeures d’échec thérapeutique. Cet article présente un protocole bioinformatique pour l’analyse transcriptomique des aspirateurs de moelle osseuse. L’objectif principal du protocole est de fournir un flux de travail standardisé pour identifier les signatures moléculaires associées à la progression de la maladie et à la résistance au traitement dans la LMA rechutée. La chaîne détaille les procédures informatiques pour comparer des échantillons de moelle osseuse non appariés, démontrées à partir de données de séquençage de cinq cas nouvellement diagnostiqués et de quatre cas rechutés. Cette méthode décrit les étapes essentielles pour traiter les données de séquençage de l’ARN, réaliser l’analyse de l’expression génique différentielle et réaliser des évaluations fonctionnelles en aval. L’application de ce flux de travail a permis d’identifier 2 025 gènes différemment exprimés (DEG), dont FOXC1, HOXA11, HOXA11-AS et AXL, comme transcripts candidats associés à la rechute dans ce jeu de données représentatif. Les analyses fonctionnelles et en réseau ont également priorisé les ensembles de gènes et les hubs d’interaction liés à la signalisation de petites GTPase, à la signalisation inflammatoire, aux interactions de matrices extracellulaires et aux processus biosynthétiques de l’ARN. Dans l’ensemble, cette méthodologie fournit un pipeline computationnel reproductible pour cartographier les signatures transcriptomiques associées à l’AML rechuté et pour générer des hypothèses nécessitant une validation expérimentale ultérieure.

Introduction

La leucémie myéloïde aiguë (LMA) est un groupe de néoplasies malignes clonales provenant de cellules souches/progénitrices hématopoïétiques, caractérisées par une prolifération anormale de cellules myéloïdes immatures dans la moelle osseuse et la suppression de la différenciation hématopoïétique ^1,2,3. Bien que la chimiothérapie d’induction standard actuelle (comme la cytarabine combinée aux anthracyclines) puisse induire une rémission complète chez la plupart des patients, le taux de rechute reste jusqu’à 50 % à 70 %, et le pronostic pour les patients en rechute est significativement compromis ^4,5,6. La résistance acquise à la chimiothérapie est associée à l’échec du traitement de la LMA, nécessitant une analyse approfondie des signatures moléculaires liées à ce processus afin d’améliorer les stratégies thérapeutiques et les taux de survie des^{patients 7,8}.

Dans l’ensemble de la littérature, des études existantes indiquent que la résistance à la chimiothérapie dans la LMA ne se limite pas à la hausse des pompes d’efflux médicamenteux ou à un métabolisme anormal des médicaments, mais est également liée au maintien de la survie des cellules souches leucémiques (CSL), à la formation de phénotypes de transition épithéliale-mésenchymatose (EMT) dans la niche hématologique, et aux interactions avec le microenvironnement de la moelle^{osseuse 9. 10,11}. Par exemple, les populations de LSC présentent une forte capacité d’auto-renouvellement et une forte quiescence, ce qui est associé à une résistance inhérente aux agents chimiothérapeutiques spécifiques au cycle^{cellulaire 12}. De plus, la régulation à la hausse des récepteurs tyrosine kinases, telles que l’AXL, a été associée à une résistance dans l’AML FLT3-ITD+, survenant parallèlement à l’activation des voies PI3K/AKT et MAPK et à des capacités anti-apoptotiques^{accrues 13,14}. La reprogrammation métabolique et le remodelage épigénétique ont également été reconnus comme des axes régulateurs importants dans la formation de résistance. Les preuves suggèrent que les cellules LMA lors de la rechute pourraient s’adapter au stress oxydatif induit par la chimiothérapie et aux dommages de l’ADN grâce à une activité accrue de phosphorylation oxydative (OXPHOS), à la modulation des rapports NAD⁺/NADH et à des états modifiés de modification^{des histones} 15,16,17. Des facteurs inflammatoires dans le microenvironnement de la moelle osseuse, tels que l’IL-6 et le CXCL8, sont également associés à la survie des LSC et à la résistance à la chimiothérapie, survenant souvent en coordination avec l’activation des voies de signalisation STAT3/^{NF-κB 11,18}.

Malgré ces mécanismes reconnus, les modifications transcriptomiques associées à la rechute de la LMA et à la chimiorésistance restent incomplètement caractérisées, notamment lorsque les signaux métaboliques, épigénétiques et liés au microenvironnement de la moelle osseuse sont évalués ensemble dans des échantillons cliniques. Ce flux de travail répond à ce besoin en priorisant les DEG candidats, les voies et les réseaux réglementaires associés à la transition du diagnostic initial à la rechute clinique. Cette approche intègre l’analyse d’expression différentielle, l’analyse d’enrichissement des ensembles de gènes (GSEA) et la construction de réseaux d’interaction protéine-protéine (PPI) pour cartographier la reprogrammation transcriptionnelle à l’échelle du système et pour nommer des candidats pour une validation mécanistique ultérieure.

L’objectif global de cette méthode est de présenter un pipeline bioinformatique standardisé et reproductible pour comparer les transcriptomes en vrac des échantillons de moelle osseuse de LMA nouvellement diagnostiqués et ceux ayant rechuté. La raison d’être de cette technique in silico est sa capacité à capturer des événements transcriptionnels non biaisés à l’échelle du génome, dépassant les limites des analyses à voie unique pour prioriser des réseaux régulateurs complexes et multidimensionnels. Cette technique offre des avantages significatifs par rapport aux méthodes alternatives, telles que les microarrays ou les panels multiplex qPCR ciblés, en offrant une plage dynamique plus élevée, la capacité de détecter de nouveaux transcrits, et une quantification précise de l’expression génique sans les limitations des sondes^{préconçues 19,20}. Pour déterminer si cette méthode est appropriée à leur application, les lecteurs doivent noter que ce pipeline est spécifiquement conçu pour les chercheurs traitant des données de séquençage d’ARN en vrac provenant de cohortes cliniques appariées ou non appariées, telles que les aspirations tissulaires. Il convient à l’identification de signatures larges associées à la résistance et des réseaux de régulation candidats, tandis que les chercheurs nécessitant une résolution spécifique au type cellulaire ou spatiale devraient utiliser des flux de travail complémentaires de séquençage unicellulaire ou spatial. En fin de compte, ce protocole informatique permet la priorisation des gènes candidats, des voies et des réseaux régulateurs pour des investigations expérimentales ultérieures.

Protocole

Toutes les méthodes impliquant l’utilisation d’échantillonnages de tissus humains ont été réalisées conformément aux directives institutionnelles et à la Déclaration d’Helsinki (révisée en 2013). Les échantillons cliniques de moelle osseuse ont été obtenus avec l’approbation du comité d’éthique institutionnel (Agréments n° TY-ZKY2024-116-01 et TY-ZKY2024-116-02).

1. Collecte clinique des échantillons et classification des patients

1. Sélectionnez des échantillons d’aspiration de moelle osseuse provenant de patients formellement diagnostiqués avec une leucémie myéloïde aiguë (LMA) selon les critères de classification de l’Organisation mondiale de la santé (OMS).
2. Appliquer des critères spécifiques d’inclusion et d’exclusion lors de la sélection des patients afin d’assurer l’homogénéité et la reproductibilité des cohortes. Inclure les patients adultes atteints de LMA primaire, et exclure les patients atteints de LMA secondaire, leucémie promyélocytaire aiguë ou antécédents d’autres cancers (Tableau 1).
3. Attribuer l’échantillon collecté au groupe nouvellement diagnostiqué si le patient présente une LMA non traitée au moment du diagnostic clinique initial.
4. Attribuer l’échantillon prélevé au groupe rechuté si le patient présente une réapparition de blasts leucémiques dans le sang périphérique ou de blasts supérieurs à 5 % dans la moelle osseuse après une rémission complète documentée.
5. Prélever les échantillons de moelle osseuse résiduelle non identifiée immédiatement après la procédure clinique d’aspiration de moelle osseuse de routine.
   REMARQUE : Dans la démonstration de ce protocole spécifique, neuf échantillons consécutifs (cinq nouvellement diagnostiqués et quatre rechutés) ont été prélevés entre septembre 2024 et septembre 2025. Comme cela n’utilisait que des échantillons cliniques résiduels désidentifiés, le comité d’éthique a levé l’exigence d’un consentement éclairé écrit.
6. Traiter immédiatement l’aspiration de moelle osseuse collectée pour la préservation de l’ARN en utilisant une méthode standard de lyse phénol-guanidinium²¹.
   1. Transférez l’aspirate fraîchement prélevée de la moelle osseuse dans un tube de collecte contenant un anticoagulant. Secouez le tube pour bien mélanger l’aspirat et l’anticoagulant.
   2. Extraire un volume mesuré de l’échantillon anticoagulé et l’ajouter à un réactif commercial de lyse phénol-guanidinium. Maintenez un rapport volumique de 3 parties de réactif de lyse pour 1 partie d’échantillon.
      ATTENTION : Le réactif de lyse phénol-guanidinium contient des produits chimiques toxiques et corrosifs pouvant provoquer de graves brûlures et des dommages tissulaires. Effectuez toute la manipulation des réactifs à l’intérieur d’une hotte chimique tout en portant un équipement de protection individuelle approprié.
   3. Secouez vigoureusement le tube pour homogénéiser complètement l’échantillon et le réactif de lyse. Assurez-vous que le mélange est bien mélangé et vérifiez qu’aucun caillot visible ne reste dans la solution.
   4. Congelez immédiatement l’échantillon homogénéisé en immergeant le tube dans de l’azote liquide.
      ATTENTION : L’azote liquide est extrêmement froid et peut provoquer de graves engelures au contact. Portez des gants cryogéniques et une visière complète lors de la manipulation de l’azote liquide.
7. Transférer les échantillons congelés instantanément dans un congélateur à -80 °C pour un stockage à long terme avant la chaîne d’isolement de l’ARN en aval et du séquençage du transcriptome. Cela représente un point sûr où l’expérience peut être mise en pause puis relancée ultérieurement.
   REMARQUE : Le flux de travail présenté dans ce protocole se concentre entièrement sur la génération de signatures de résistance computationnelle. Aucune validation expérimentale indépendante, telle que la PCR quantitative en temps réel (RT-qPCR), n’a été réalisée sur les principaux gènes différenciellement exprimés identifiés via ce pipeline spécifique.

2. Contrôle de la qualité de l’ARN et préparation des bibliothèques

1. Évaluez l’intégrité de l’ARN à l’aide d’un système d’électrophorèse capillaire microfluidique. Pour ce flux de travail représentatif, incluez des échantillons d’ARN avec un nombre d’intégrité ARN (RIN) ≥ 6,0, un rapport A260/280 entre 1,8 et 2,1, et aucun pic de dégradation visible. Enregistrez les ratios mesurés de RIN et de pureté pour chaque échantillon avant la préparation à la bibliothèque.
2. Apportez 1 μg d’ARN total par échantillon pour la préparation de la bibliothèque. Purifiez l’ARNm à partir de l’ARN total en utilisant des billes magnétiques liliées à l’oligo poly-T pour enrichir les transcrits à queue polyA.
3. Fragmentez l’ARNm enrichi à l’aide de cations divalents. Incubez le mélange à 94 °C pendant 15 minutes dans un tampon de réaction de synthèse de première brin 5X.
4. Synthétiser l’ADNc de première brin à l’aide d’amorces hexamères aléatoires et d’une transcriptase inverse dépourvue d’activité RNase H.
5. Dégradez le brin modèle d’ARN à l’aide de la RNase H. Synthétissez le cDNA du second brin en utilisant la polymérase ADN I et les dNTP dans un système de réaction de 20 μL.
6. Incubez la réaction de synthèse du second brin à 16 °C pendant 1 heure. Centrifugez brièvement le mélange de réaction à 2 000 x g pour recueillir le liquide au fond du tube.
7. Convertir les débordements restants en extrémités émoussées en utilisant des activités d’exonucléase et de polymérase. Adénylate les extrémités 3' des fragments d’ADN et ligate les adaptateurs avec des structures en boucle en épingle pour préparer l’hybridation.
8. Purifiez les fragments de la bibliothèque à l’aide de billes d’immobilisation magnétiques en phase solide réversible pour sélectionner préférentiellement des fragments d’ADNc de 370 à 420 pb de long.
9. Effectuez des lavages à l’éthanol pendant la purification des billes. Centrifugez les tubes à 2 000 x g pendant 30 s pour collecter et éliminer complètement tout éthanol résiduel avant l’élution finale.
10. Effectuer une amplification PCR à l’aide d’une ADN polymérase haute fidélité, des amorces PCR universelles et des amorces à index spécifiques à chaque échantillon.
11. Exécutez le profil thermique PCR avec une dénaturation initiale à 98 °C pendant 30 s. Suivez 12 cycles à 98 °C pendant 10 s, 60 °C pendant 30 s et 72 °C pendant 30 s, terminant par une extension finale à 72 °C pendant 5 minutes.
12. Purifiez à nouveau les produits PCR avec les billes magnétiques. Appliquez les mêmes paramètres de centrifugation que l’étape 2.9 pour obtenir la bibliothèque finale.
13. Quantifiez la concentration initiale de la bibliothèque à l’aide d’un fluoromètre. Diluez la bibliothèque finale à une concentration de 1,5 ng/μL.
14. Mélangez soigneusement la bibliothèque diluée. Centrifugez le mélange à 10 000 x g pendant 1 minute à 4 °C pour éliminer les débris restants avant l’analyse finale.
15. Évaluez la taille de l’insert à l’aide du système d’électrophorèse capillaire microfluidique.
16. Quantifiez avec précision la concentration effective de la bibliothèque via la PCR quantitative en temps réel (qRT-PCR) après avoir confirmé que la taille de l’insert répond aux attentes. Assurez-vous que la concentration dépasse 1,5 nM afin de garantir la stabilité de la bibliothèque et la qualité du séquençage.
    REMARQUE : Cela représente un point sûr où l’expérience peut être mise en pause. Les bibliothèques préparées peuvent être stockées à -20 °C jusqu’au clustering et au séquençage.

3. Regroupement et séquençage du transcriptome

1. Effectuez le regroupement des échantillons codés par indice sur un système automatisé de génération de cluster. Utilisez un kit commercial de cluster à extrémités jumelées selon les instructions du fabricant.
2. Séquencez les préparations de la bibliothèque sur une plateforme de séquençage à haut débit après la génération réussie du cluster. Générez 150 lectures par paires de bases (pb) à l’extrémité de la b.

4. Contrôle de la qualité des données et cartographie des lectures

1. Évaluez la qualité des données brutes (format FASTQ) en utilisant fastp v0.23.2 pour le contrôle qualité et le filtrage en lecture brute. Enregistrez les paramètres de la ligne de commande dans un journal d’analyse. Dans ce flux de travail, des lectures propres étaient générées en supprimant les lectures contenant des adaptateurs, les lectures contenant des séquences poly-N, et les lectures de faible qualité en utilisant des réglages de filtrage identiques entre les échantillons. Une commande représentative en binôme est fournie dans le Fichier Supplémentaire 1.
2. Traitez les lectures brutes via un logiciel de prétraitement automatisé. Obtenez des lectures claires en supprimant les lectures contenant des adaptateurs, celles contenant des séquences poly-N et les lectures de faible qualité. Utilisez des paramètres de filtrage identiques pour tous les échantillons et enregistrez le numéro de lecture conservé, Q20, Q30 et le contenu GC après filtrage.
3. Calculez le contenu Q20, Q30 et GC des données propres. Définir les variables potentielles du lot avant l’analyse en aval, y compris la date de collecte de l’échantillon, la date d’extraction de l’ARN, le lot de préparation de la bibliothèque, la voie de séquençage et la course de séquençage.
4. Évaluer les effets des lots par ACP et analyse de corrélation échantillon à échantillon en utilisant des valeurs d’expression normalisées. Si les échantillons se regroupent principalement par variables techniques plutôt que par état clinique, documentez la variable affectée et incluez-la comme covariée dans la formule de conception d’expression différentielle ou appliquez une méthode établie d’ajustement par lots avant la visualisation en aval.
5. Acquérez le génome de référence (Homo sapiens, GRCh38) et les fichiers d’annotation génétique Ensembl correspondants pour libérer 109 pour l’alignement de lecture.
6. Construis l’index du génome de référence en utilisant HISAT2 v2.0.5.
7. Alignez les lectures claires à l’extrémité jumelée au génome de référence en utilisant HISAT2 v2.0.5. Utilisez cette approche d’alignement conscient des épissures pour générer une base de données de jonctions de splince basée sur le fichier d’annotation du modèle génétique.

5. Prédiction de nouvelles transcriptions et quantification de l’expression génique

1. Assemblez les lectures mappées de chaque échantillon à l’aide de StringTie v1.3.3b selon une approche basée sur la référence. Utilisez cet outil pour assembler et quantifier des transcrits de longueur complète représentant plusieurs variantes de sépissage pour chaque locus génique.
2. Comptez le nombre de lectures mappées à chaque gène en utilisant featureCounts v1.5.0-p3. Utilisez la matrice brute de décompte des entiers comme entrée pour l’analyse d’expression différentielle en aval.
3. Configurez featureCounts v1.5.0-p3 pour les données de séquençage des extrémités appariées en utilisant l’option des extrémités appariées (par exemple, -p). Fournissez le fichier d’annotation GTF GRCh38 téléchargé pour définir les limites correctes des caractéristiques génomiques.
4. Calculez les fragments par kilobase de transcrit par million de lectures cartographiées (FPKM) pour chaque gène. Utilisez uniquement les valeurs FPKM pour la visualisation descriptive, l’ACP, l’affichage de la carte thermique et les résumés d’expressions exploratoires ; n’utilisent pas les valeurs FPKM comme matrice d’entrée pour les tests d’expression différentielle DESeq2.

6. Analyse différentielle de l’expression génique

1. Effectuer une analyse d’expression différentielle entre les groupes nouvellement diagnostiqués et les groupes rechutés en utilisant R v3.5.0 et le package DESeq2 R v1.20.0. Importez la matrice brute de décompte de lectures générée à l’étape 5.2 dans l’environnement R, et conservez les valeurs FPKM uniquement pour la visualisation et les analyses exploratoires.
2. Construisez l’objet de jeu de données spécialisé requis par le paquet d’analyse. Exécutez la commande spécifique (par exemple, DESeqDataSetFromMatrix()) pour lier la matrice de données de comptage avec la table de métadonnées d’exemple correspondante.
3. Définissez la formule de conception expérimentale au sein de l’objet logiciel. Spécifiez l’état clinique (nouvellement diagnostiqué versus rechute) comme variable principale de comparaison (par exemple, conception = ~ condition). Si une variable technique du lot est identifiée à l’étape 4.3 et n’est pas complètement confondue avec l’état clinique, incluez-la dans la formule de conception (par exemple, conception = ~ lot + condition).
4. Exécutez la fonction d’analyse différentielle d’expression centrale (par exemple, DESeq()). Permettre au logiciel d’effectuer automatiquement l’estimation du facteur de taille, l’estimation de dispersion et le test de test binomial^{Wald 22} négatif.
5. Extraire la table des résultats en utilisant la fonction d’extraction des résultats (par exemple, résultats()). Spécifiez l’argument du contraste pour définir la comparaison exacte (rechute versus nouveau diagnostiqué).
6. Ajustez les valeurs p résultantes pour contrôler le taux de fausses découvertes. Utilisez la procédure intégrée Benjamini et Hochberg appliquée automatiquement par le logiciel²³.
7. Filtrez la table des résultats extraits pour isoler les gènes différentiellement exprimés (DEG) significatifs. Attribuer à tout gène avec une valeur P ajustée < 0,05 et un logarithmide absolu₂ fois changement > 1 comme exprimé de manière significativement différentielle.

7. Analyse d’enrichissement par ontologie génique (GO)

1. Effectuer une analyse d’enrichissement par ontologie génétique (GO) des DEG identifiés à l’aide de clusterProfiler v3.8.1 et org. Hs.eg.db v3.6.0. Entrez la liste des identifiants de gènes Entrez correspondant aux DEG significatifs identifiés à l’étape 6.7.
2. Exécutez la fonction d’enrichissement GO (par exemple, enrichGO()). Spécifiez les paramètres requis, y compris la base de données d’organismes de fond appropriée (par exemple, OrgDb = org. Hs.eg.db), le domaine ontologique spécifique (Processus biologique, composante cellulaire ou fonction moléculaire), et le seuil ajusté de la valeur P (0,05).
3. Assurez-vous que l’algorithme applique les corrections nécessaires lors du calcul de l’enrichissement. Confirmez que le logiciel corrige en interne le biais de longueur du gène et ajuste les valeurs P en utilisant la méthode²⁴ de Benjamini et Hochberg.
4. Considérons les termes GO avec une valeur P corrigée inférieure à 0,05 comme significativement enrichis. Générez un diagramme de points ou un graphique à barres en utilisant les fonctions de visualisation intégrées du package pour afficher les termes GO enrichis les plus hauts.

8. Encyclopédie Kyoto des gènes et génomes (KEGG) Analyse de l’enrichissement des voies

1. Utilisez une ressource de base de données complète dédiée à la compréhension des fonctions de haut niveau du système biologique afin d’identifier les voies dérégulées. Préparez la même liste des identifiants DEG Entrez principaux utilisés à l’étape 7.1.
2. Exécutez la fonction d’enrichissement KEGG (par exemple, enrichKEGG()) dans le paquet d’annotation fonctionnelle R.
3. Définissez les paramètres critiques dans l’appel de fonction. Définir strictement le code de l’organisme sur humain (par exemple, organisme = « possède ») et définir la méthode d’ajustement de la valeur P (par exemple, méthode pAdadjust = « BH »).
4. Extraire les voies KEGG statistiquement significatives. Filtrez la sortie pour ne conserver que les voies montrant une valeur P corrigée inférieure à 0,05.
5. Visualisez les voies KEGG les plus riches. Utilisez les fonctions de tracé intégrées (par exemple, dotplot()) pour cartographier la signification statistique et le nombre de gènes associés à chaque voie.

9. Analyse d’enrichissement des ensembles géniques (GSEA)

1. Préparez la liste génétique préclassée requise pour l’analyse. Calculez la métrique de classement pour tous les gènes exprimés en utilisant le signe -log10 (valeur P) multiplié par le signe du changement de_{log 2} multiplié par l’analyse différentielle d’expression.
2. Lancez une installation locale du logiciel GSEA v4.2.3 du Broad Institute. Saisissez la nouvelle liste génétique pré-classée dans l’interface^{logicielle 25}.
3. Téléchargez les ensembles de gènes prédéfinis requis. Acquérir les ensembles de données Gene Ontology (GO) et Kyoto Encyclopedia of Genes and Genomes (KEGG) à partir de la base de données Molecular Signatures (MSigDB, version 7.5.1)²⁶.
4. Configurez les paramètres logiciels pour effectuer le test d’enrichissement statistique. Réglez le nombre de permutations à 1 000 et sélectionnez le type de permutation comme « gene_set ».
5. Exécutez l’algorithme d’analyse pour déterminer si les ensembles de gènes prédéfinis présentent une différence statistiquement significative et concordante entre les états biologiques nouvellement diagnostiqués et ceux qui ont rechuté.
6. Évaluer la signification statistique des profils d’enrichissement générés. Définir des ensembles de gènes significatifs en utilisant des seuils stricts : une valeur absolue du score d’enrichissement normalisé (NES) > 1,0, une valeur P nominale < 0,05, et une valeur q du taux de fausse découverte (FDR) < 0,25.

10. Analyse du réseau d’interaction protéine-protéine (IPP)

1. Accédez à la base de données STRING pour les interactions protéine-protéine connues et prédites. Dans ce flux de travail, l’analyse PPI a été réalisée avec STRING v11.5²⁷.
2. Entrez la liste des identifiants de gènes Entrez ou des symboles officiels des gènes significatifs exprimés différemment (identifiés à l’étape 6.7) dans l’interface de recherche de la base de données. Sélectionnez Homo sapiens comme organisme cible.
3. Configurez les paramètres de construction du réseau pour garantir que des interactions de haute qualité soient récupérées. Fixez le score minimum requis d’interaction à un seuil de confiance élevé (score > 0,700).
4. Exportez les données d’interaction du réseau résultantes dans un annuaire local. Enregistrez la carte d’interaction sous forme de fichier tabulaire standard (par exemple, au format TSV).
5. Importez les données d’interaction exportées dans Cytoscape v3.9.1 pour la visualisation et l’analyse du^{réseau 28}.
6. Filtrez le réseau construit pour améliorer la clarté de la visualisation et mettre en évidence les principaux pôles réglementaires. Retirez tout nœud déconnecté ou gène orphelin qui ne présente pas d’interactions continues atteignant le seuil de confiance établi.

Résultats

Validation de la cohorte clinique et du séquençage

La réussite de l’exécution du protocole d’extraction et de préparation de bibliothèque d’ARN en amont (Figure 1) a été confirmée par le séquençage des métriques de rendement et de qualité. Dans ce jeu de données représentatif, des échantillons de moelle osseuse provenant de cinq patients nouvellement diagnostiqués avec LMA et de quatre patients atteints de LMA rechuté ont fourni en moyenne environ 6,0 Go de données brutes par échantillon. L’évaluation du contrôle qualité (Tableau 2) a confirmé que la qualité de base et la profondeur de lecture respectaient les seuils requis pour l’analyse bioinformatique en aval9. Une faible intégrité de l’ARN (par exemple, RIN < 6,0), de faibles taux de cartographie ou un biais élevé de dégradation des transcrits représenteraient une qualité d’entrée sous-optimale et compromettraient la fiabilité de l’analyse différentielle d’expression en aval.

Variance transcriptomique globale et ACP

Pour évaluer la variance transcriptomique globale et inspecter le regroupement clinique, l’ACP a été réalisée sur les données d’expression normalisée. Dans cet ensemble de données représentatif, les groupes nouvellement diagnostiqués et rechutes ont montré une séparation dans l’espace bidimensionnel (Figure 2A)²⁰, PC1 et PC2 représentant respectivement 23,82 % et 18,75 % de la variance totale. Les diagrammes de Venn de la Figure 2B,C fournissent un résumé descriptif supplémentaire des gènes détectés parmi les échantillons des groupes nouvellement diagnostiqués et rechutés, soutenant les vérifications de reproductibilité au niveau de l’échantillon avant l’analyse d’expression différentielle en aval. Parce que la cohorte était petite et non appariée, la séparation des ACP a été interprétée comme un résultat illustratif du flux de travail plutôt que comme une preuve définitive d’une biologie spécifique à l’état de la maladie.

Analyse du gène d’expression différentielle (DEG)

L’application des seuils de protocole établis (|_{log 2}FC| ≥ 1 et la valeur P ajustée ≤ 0,05) à la sortie DESeq2 a permis d’identifier 2 025 DEGs, comprenant 772 gènes à la hausse et 1 253 à la baisse dans le groupe rechuté (Figure 3A). Les relevés de notes candidats avec une forte variation comprenaient FOXC1 (log₂FC = 7,55, P = 4,92 x ^10-5), HOXA11 (log₂FC = 7,76), HOXA11-AS (log₂FC = 7,23) et AXL (log₂FC = 3,50), ainsi que RHOB, PTX3 et CXCL8 régulés à la baisse. La littérature existante relie plusieurs de ces gènes à la stempillage, à la signalisation ou à la réponse thérapeutique de la^{LMA 13,29} ; cependant, le flux de travail actuel ne les identifie que comme des relevés de notes de candidats associés à une rechute. Tout rôle mécaniste définitif dans la résistance clinique nécessite une validation fonctionnelle indépendante ultérieure.

Enrichissement fonctionnel et des voies (GO, KEGG et GSEA)

Le protocole d’annotation fonctionnelle a mappé les DEG vers des systèmes biologiques plus larges. L’analyse GO a identifié un enrichissement de termes liés à la transduction de signaux médiée par la GTPase petite, au transport d’ions métalliques et à l’assemblage de la chromatine (Figure 4A–C). La cartographie des voies KEGG a identifié des associations avec les interactions récepteur ECM et cytokine-récepteur (Figure 4D). GSEA a montré un enrichissement des processus biosynthétiques de l’ARN dans le groupe rechuté et un enrichissement des voies du métabolisme énergétique dans le groupe nouvellement diagnostiqué (Figure 5A). Ces résultats d’enrichissement fournissent une feuille de route descriptive des ensembles génétiques modifiés et doivent être interprétés comme des associations génératrices d’hypothèses plutôt que comme des facteurs prouvés de rechute.

Construction de réseaux d’interaction protéine-protéine (IPP)

Le réseau STRING initial contenait 56 nœuds et 193 interactions. Après suppression des nœuds déconnectés ou orphelins, le sous-réseau Cytoscape affiché contenait 42 nœuds et 136 interactions (Figure 5B). L’analyse modulaire réseau a donné la priorité à TP53, CCL2, CXCL8 et IL6 comme centres mathématiques centraux avec le plus grand nombre d’interactions. Parce que le réseau PPI repose sur des scores d’interaction prédits dans la base de données (par exemple, score ATF3 : 0,982), l’identification du hub doit être interprétée comme une priorisation de cible pour de futures études empiriques plutôt que comme une preuve directe d’une évasion d’apoptose médiée par p53 ou d’autres mécanismes de résistance.

Les données brutes de séquençage ARN générées dans ce protocole ont été déposées dans le dépôt Figshare et sont accessibles au public via le DOI suivant : https://doi.org/10.6084/m9.figshare.30655814. Les données traitées et les fichiers d’analyse associés sont inclus dans l’article et ses documents complémentaires. Les paramètres représentatifs en ligne de commande et les paramètres d’analyse utilisés pour reproduire le flux de travail computationnel sont fournis sous forme de fichier supplémentaire 1. Toutes les données soutenant les résultats de cette étude sont disponibles sans restriction.

Identification du patient	Âge (Années)	Sexe	Mutations moléculaires	Survie/Suivi (Mois)	Statut clinique
R_AML_1	70	Homme	FLT3-ITD (+)	22	Décédé
R_AML_2	29	Féminin	NPM1 (+)	11	Vivant
R_AML_3	40	Homme	CEBPA (+)	17	Vivant
R_AML_4	55	Féminin	Triple négatif*	24	Décédé

Tableau 1 : Caractéristiques démographiques et cliniques des patients du groupe LMA rechuté (R_AML). Le tableau 1 résume les caractéristiques démographiques et cliniques de la cohorte de LMA rechute utilisée dans l’analyse représentative, y compris les caractéristiques cliniques au niveau du patient pertinentes pour l’interprétation du flux de travail transcriptomique.

Exemple	Bibliothèque	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

Tableau 2 : Résumé de la qualité des données. Le tableau 2 rapporte les métriques de qualité du séquençage pour chaque échantillon, incluant le rendement en lecture, la qualité de base, le contenu GC et les informations de contrôle qualité liées à la cartographie utilisées pour déterminer si les échantillons étaient adaptés à une analyse en aval.

Figure 1 : Flux de travail du protocole. Le flux de travail résume les principales étapes expérimentales et computationnelles, incluant la collecte d’échantillons cliniques, le contrôle qualité ARN, la préparation et le séquençage des bibliothèques, le traitement et l’alignement de la lecture, la quantification des transcriptions, l’analyse d’expression différentielle, l’enrichissement GO/KEGG, la GSEA et la construction du réseau PPI. Veuillez cliquer ici pour voir une version agrandie de cette figurine.

Figure 2 : Analyse quantitative des échantillons. (A) L’analyse en composantes principales (ACP) a été réalisée pour évaluer les différences entre groupes et la reproductibilité de l’échantillon au sein du groupe. L’ACP a été réalisée à l’aide de méthodes algébriques linéaires basées sur des valeurs d’expression génique normalisées pour tous les échantillons. (B, C) Diagrammes de Venn montrant les gènes détectés à travers les échantillons des groupes AML et R_AML, respectivement. Les régions restreintes par échantillon indiquent les gènes détectés dans des échantillons individuels, tandis que les zones qui se chevauchent représentent des gènes couramment détectés sur deux échantillons ou plus. Veuillez cliquer ici pour voir une version agrandie de cette figurine.

Figure 3 : Analyse différentielle de l’expression génique. (A) Graphique en barres montrant le nombre de gènes (DEG) exprimés différemment entre groupes de comparaison, identifiés par DESeq2 avec des seuils de valeur P ajustée ≤ 0,05 et |_{log 2}FoldChange| ≥ 1. (B) Parcelle volcanique des DEG. L’axe des x représente les valeurs de_{FoldChange de log 2}, et l’axe des y représente -log10 (valeur P). Les lignes bleues pointillées indiquent les seuils utilisés pour la sélection des DEG. (C) Carte thermique de regroupement hiérarchique des DEG. L’axe des x indique les noms des échantillons, et l’axe des y montre les valeurs d’expression normalisées des DEG. Veuillez cliquer ici pour voir une version agrandie de cette figurine.

Figure 4 : Analyse d’enrichissement fonctionnel des gènes exprimés différemment. (A) Parcelle du bar d’enrichissement GO. L’axe des x représente les termes GO, et l’axe des y montre la signification d’enrichissement, exprimée en -log10(padj). Les couleurs représentent BP (processus biologique), CC (composante cellulaire) et MF (fonction moléculaire). (B) Graphique de bulle d’enrichissement GO. L’axe des x représente le rapport des DEG annotés à chaque terme GO par rapport au nombre total de DEG, et l’axe des y indique les termes GO. La taille des bulles correspond au nombre de gènes annotés, et les gradients de couleur représentent la signification de l’enrichissement. (C) Parcelle du bar d’enrichissement KEGG. L’axe des x représente les voies KEGG, et l’axe des y indique la signification de l’enrichissement. (D) Complot de la bulle d’enrichissement de KEGG. La taille de la bulle indique le nombre de gènes annotés, et les gradients de couleur reflètent la signification de l’enrichissement. Veuillez cliquer ici pour voir une version agrandie de cette figurine.

Figure 5 : Enrichissement GSEA et analyse du réseau d’interaction protéine-protéine (PPI). (A) Graphique en barres montrant les scores d’enrichissement normalisés (NES) pour certains ensembles génétiques significatifs. Les valeurs positives de NES indiquent un enrichissement dans le groupe R_AML, tandis que des valeurs négatives de NES indiquent un enrichissement dans le groupe nouvellement diagnostiqué avec LMA. (B) Réseau d’interaction protéine-protéine (IPP). Chaque nœud représente une protéine, et chaque arête désigne une interaction entre des protéines connectées. Veuillez cliquer ici pour voir une version agrandie de cette figurine.

Discussion

Étapes critiques du protocole

La réussite de ce flux de travail bioinformatique repose sur plusieurs étapes cruciales. Premièrement, la congélation immédiate et la lyse correcte de l’aspirate de la moelle osseuse (Étape 1.6) sont primordiales, car le microenvironnement de la moelle osseuse est riche en ribonucléases capables de dégrader rapidement l’intégrité^{transcriptomique 30}. Pendant la phase de calcul, la configuration correcte de la formule de conception expérimentale dans le package DESeq2 (étape 6.3) est cruciale pour une expression différentielle précise, en particulier lorsque l’état clinique (nouvellement diagnostiqué versus rechute) est comparé tout en considérant les variables de confusion potentielles. Enfin, appliquer des seuils rigoureux de taux de fausse découverte (FDR) lors de l’analyse d’enrichissement des ensembles géniques (GSEA) (étape 9.6) constitue un point de contrôle statistique critique pour éviter la surinterprétation des réseaux fonctionnels faussement positifs.

Modifications et dépannage

Un défi courant dans cette méthode est la présence d’effets batch, qui surviennent fréquemment lors de la collecte et du séquencage d’échantillons cliniques sur des périodes prolongées. Les variables de lot doivent être définies avant analyse, y compris la date de collecte de l’échantillon, la date d’extraction de l’ARN, le lot de préparation de la bibliothèque, la voie de séquençage et la course de séquençage. Si l’analyse PCA ou corrélation d’échantillon révèle un regroupement basé sur la date de séquençage ou une autre variable technique plutôt que sur le phénotype clinique, les utilisateurs doivent modifier le protocole en incluant la variable batch dans la formule de conception d’expression différentielle lorsque cela est statistiquement réalisable ou en appliquant des algorithmes de correction par lots, tels que ComBat ou SVA, avant la visualisation³¹. Si ce protocole est appliqué au sang total plutôt qu’aux aspirateurs de moelle osseuse, une modification essentielle est l’inclusion d’une étape de déplétion de l’ARNm de la globine lors de la préparation de la bibliothèque afin d’empêcher les transcrits de globine très abondants de monopoliser la profondeur de lecture du séquençage. Les versions logicielles et les paramètres principaux du flux de travail représentatif ont été complétés comme suit : 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, séquençage en binôme à 150 pb, GSEA v4.2.3 avec 1 000 permutations de jeux de gènes, MSigDB v7.5.1, STRING v11.5 avec interactions à haute confiance, et Cytoscape v3.9.1. Les paramètres représentatifs de la ligne de commande et les paramètres d’analyse sont fournis dans le Fichier Supplémentaire 1.

Limites de la méthode

Bien que global, ce protocole présente des limites méthodologiques inhérentes. Premièrement, il utilise le séquençage d’ARN en masse, qui capture le profil transcriptomique moyen de l’ensemble de l’aspirate de la moelle osseuse et ne possède pas de résolution spatiale unicellulaire. Par conséquent, le flux de travail ne peut pas déterminer si une signature associée à la rechute surélevée provient des cellules souches leucémiques, des cellules stromales, des cellules immunitaires ou des changements de composition cellulaire³². Deuxièmement, le jeu de données représentatif est petit (n = 9) et non apparié, ce qui limite la robustesse statistique et empêche une inférence causale définitive. Troisièmement, le flux de travail est purement in silico. Il génère des centres régulateurs et des voies de signalisation candidats, mais il ne peut pas valider indépendamment leur nécessité fonctionnelle en chimiorésistance sans validation expérimentale in vitro ou in vivo orthogonale.

Des études récentes sur la génomique unicellulaire et la cellule unique ont élargi le cadre de référence de la LMA en résolvant l’hétérogénéité des états cellulaires, l’architecture clonale et l’évolution associée à la thérapie à une résolution plus^élevée 33,34,35,36. Ces approches sont complémentaires au flux de travail de RNA-seq en vrac décrit ici : le séquençage en vrac offre une stratégie de dépistage pratique et rentable des signatures transcriptomiques au niveau de la cohorte, tandis que les méthodes monocellulaires et multi-omiques peuvent être utilisées dans des études de suivi pour assigner des signaux candidats à des populations cellulaires malignes ou microenvironnementales spécifiques.

Importance par rapport aux méthodes existantes

Malgré ces limitations, ce pipeline transcriptomique offre des avantages par rapport aux techniques diagnostiques et analytiques alternatives. Les évaluations cliniques traditionnelles de la rechute de la LMA reposent souvent sur des panels multiplex qPCR ciblés ou une cytométrie à flux standard. Bien qu’utiles pour des diagnostics rapides, ces méthodes ciblées sont limitées par des sondes prédéfinies et ne peuvent évaluer que les marqueurs de résistance^{connus 19}. En utilisant un séquençage transcriptomique non biaisé à l’échelle du génome, associé à une analyse en réseau, ce protocole peut nommer de nouveaux transcrits et associations systémiques que les méthodes ciblées existantes pourraient négliger.

Importance et applications potentielles

La méthodologie décrite dans ce protocole est pertinente pour l’hématologie translationnelle et la médecine personnalisée car elle peut prioriser les signatures transcriptomiques associées à la rechute pour des études supplémentaires. Une application potentielle en aval est la nomination d’antigènes de surface ou de voies d’évasion immunitaire qui émergent lors de la rechute. De tels candidats pourraient orienter la conception de futures études de validation et, s’ils sont confirmés expérimentalement, pourraient contribuer au développement d’immunothérapies de nouvelle génération, y compris les stratégies cellulaires CAR-T ou^{CAR-NK 37}. Ces applications translationnelles doivent être considérées comme génératrices d’hypothèses plutôt que comme des conclusions établies à partir du jeu de données actuel.

Matériaux

Liste des matériaux utilisés dans cet article

Nom	Entreprise	Numéro de catalogue	Commentaires
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