Perfilado transcriptómico y análisis bioinformático de muestras de médula ósea para identificar firmas de resistencia a la quimioterapia en leucemia mieloide aguda

Resumen

Este protocolo presenta un flujo de trabajo bioinformático estandarizado para analizar alteraciones transcriptómicas en la leucemia mieloide aguda (LMA). El objetivo es comparar muestras de médula ósea recién diagnosticadas y recaidas, y priorizar las firmas moleculares asociadas con la resistencia a la quimioterapia y la progresión de la enfermedad para investigaciones posteriores.

Resumen

La leucemia mieloide aguda (LMA) es una malignidad hematológica altamente heterogénea en la que la recaída y la quimioresistencia adquirida siguen siendo las principales causas de fracaso del tratamiento. Este artículo presenta un protocolo bioinformático para el análisis transcriptómico de aspirantes de médula ósea. El objetivo principal del protocolo es proporcionar un flujo de trabajo estandarizado para identificar firmas moleculares asociadas con la progresión de la enfermedad y la resistencia a la terapia en LMA recaído. La cadena detalla los procedimientos computacionales para comparar muestras de médula ósea no emparejadas, demostrados mediante datos de secuenciación de cinco casos recién diagnosticados y cuatro casos recaídos. Este método describe los pasos esenciales para procesar datos de secuenciación de ARN, realizar análisis diferenciales de expresión génica y realizar evaluaciones funcionales posteriores. Aplicando este flujo de trabajo se identificaron 2.025 genes expresados diferencialmente (DEGs), incluyendo FOXC1, HOXA11, HOXA11-AS y AXL, como transcripciones candidatas asociadas a la recaída en este conjunto de datos representativo. Los análisis funcionales y de redes priorizaron además conjuntos génicos y centros de interacción relacionados con la señalización de pequeñas GTPasas, la señalización inflamatoria, las interacciones de la matriz extracelular y los procesos biosintéticos de ARN. En conjunto, esta metodología proporciona una cadena computacional reproducible para mapear firmas transcriptómicas asociadas a la LMA recaidente y para generar hipótesis que requieren validación experimental posterior.

Introducción

La leucemia mieloide aguda (LMA) es un grupo de neoplasias malignas clonales originadas en células madre/progenitoras hematopoyéticas, caracterizadas por una proliferación anormal de células mieloides inmaduras en la médula ósea y la supresión de la diferenciación hematopoyética ^1,2,3. Aunque la quimioterapia de inducción estándar actual (como la citarabina combinada con antracicinas) puede inducir remisión completa en la mayoría de los pacientes, la tasa de recaída sigue siendo tan alta como entre el 50% y el 70%, y el pronóstico para pacientes recaídos está significativamente ^{comprometido 4,5,6}. La resistencia adquirida a la quimioterapia se asocia con el fracaso del tratamiento de la LMA, lo que requiere un análisis en profundidad de las firmas moleculares vinculadas a este proceso para mejorar las estrategias terapéuticas y las tasas de supervivencia de los pacientes ^7,8.

En la literatura más amplia, los estudios existentes indican que la resistencia a la quimioterapia en LMA no se limita a la regulación al alza de bombas de eflujo de fármacos o a un metabolismo anormal de fármacos, sino que también está vinculada al mantenimiento de la supervivencia de las células madre leucemias (CSL), la formación de fenotipos similares a la transición epitelial-mesenquimatosa (EMT) dentro del nicho hematológico y las interacciones con el microambiente de la médula^{ósea 9. 10,11}. Por ejemplo, las poblaciones de LSC exhiben una alta capacidad de autorrenovación y quiescencia, lo que se asocia con una resistencia inherente a agentes quimioterapéuticos específicos del ciclo^{celular 12}. Además, la regulación al alza de las tirosinas quinasas receptoras, como AXL, se ha asociado con resistencia en FLT3-ITD+ AML, que ocurre junto con la activación de las vías PI3K/AKT y MAPK y con capacidades antiapoptóticas^{mejoradas 13,14}. La reprogramación metabólica y la remodelación epigenética también han sido reconocidas como ejes reguladores importantes en la formación de resistencias. La evidencia sugiere que las células LMA durante la recaída pueden adaptarse al estrés oxidativo inducido por quimioterapia y al daño en el ADN mediante una actividad de fosforilación oxidativa aumentada (OXPHOS), la modulación de las relaciones NAD⁺/NADH y la alteración de los estados de modificación de las^histonas 15,16,17. Los factores inflamatorios en el microambiente de la médula ósea, como IL-6 y CXCL8, también están asociados con la supervivencia a los LSC y la resistencia a la quimioterapia, ocurriendo a menudo en coordinación con la activación de las vías de señalización STAT3/^{NF-κB 11,18}.

A pesar de estos mecanismos reconocidos, los cambios transcriptómicos asociados a la recaída de la LMA y la quimioresistencia permanecen incompletamente caracterizados, especialmente cuando se evalúan señales metabólicas, epigenéticas y relacionadas con el microambiente de la médula ósea en muestras clínicas. Este flujo de trabajo responde a esta necesidad priorizando los DEG candidatos, las vías y las redes regulatorias asociadas con la transición del diagnóstico inicial a la recaída clínica. El enfoque integra análisis de expresión diferencial, análisis de enriquecimiento de conjuntos génicos (GSEA) y construcción de redes de interacción proteína-proteína (PPI) para mapear la reprogramación transcripcional a nivel de sistema y nominar candidatos para una posterior validación mecanicista.

El objetivo general de este método es presentar una cadena estandarizada y reproducible de bioinformática para comparar los transcriptomas masivos de muestras de médula ósea de LMA recién diagnosticadas frente a las recaídas. La razón de su uso de esta técnica in silico es su capacidad para capturar eventos transcripcionales imparciales y de alcance genómico, superando las limitaciones de los análisis de vía única para priorizar redes reguladoras complejas y multidimensionales. Esta técnica ofrece ventajas significativas sobre métodos alternativos, como microarrays o paneles multiplex qPCR dirigidos, al proporcionar un rango dinámico más alto, la capacidad de detectar transcritos nuevos y una cuantificación precisa de la expresión génica sin las limitaciones de las sondas^{prediseñadas 19,20}. Para determinar si este método es adecuado para su aplicación, los lectores deben tener en cuenta que esta cadena está diseñada específicamente para investigadores que procesan datos de secuenciación de ARN a granel de cohortes clínicas emparejadas o no emparejadas, como aspiraciones de tejidos. Es adecuado para identificar firmas amplias asociadas a la resistencia y redes regulatorias candidatas, mientras que los investigadores que requieren resolución específica o espacial por tipo celular tendrían que utilizar flujos de secuenciación complementarios de célula única o espacial. En última instancia, este protocolo computacional permite priorizar genes candidatos, vías y redes regulatorias para investigaciones experimentales posteriores.

Protocolo

Todos los métodos que implicaban el uso de muestreo de tejido humano se realizaron en cumplimiento de las directrices institucionales y de la Declaración de Helsinki (revisada en 2013). Las muestras clínicas de médula ósea se obtuvieron con la aprobación del comité de ética institucional (Números de aprobación TY-ZKY2024-116-01 y TY-ZKY2024-116-02).

1. Recogida de muestras clínicas y clasificación de pacientes

1. Selecciona muestras de aspiración de médula ósea de pacientes diagnosticados formalmente con leucemia mieloide aguda (LMA) según criterios de clasificación de la Organización Mundial de la Salud (OMS).
2. Aplicar criterios específicos de inclusión y exclusión durante la selección del paciente para asegurar la homogeneidad y reproducibilidad de las cohortes. Incluir pacientes adultos con LMA primaria y excluir a los pacientes con LMA secundaria, leucemia promielocítica aguda o antecedentes previos de otras malignidades (Tabla 1).
3. Asignar la muestra recogida al grupo recién diagnosticado si el paciente presenta LMA no tratado previamente en el momento del diagnóstico clínico inicial.
4. Asignar la muestra recogida al grupo recaente si el paciente muestra una reaparición de blastos leucémicos en la sangre periférica o blastos superiores al 5% en la médula ósea tras una remisión completa documentada.
5. Recoger las muestras residuales de médula ósea identificadas inmediatamente después del procedimiento clínico rutinario de aspiración de médula ósea.
   NOTA: En la demostración de este protocolo específico, se recogieron nueve muestras consecutivas (cinco recién diagnosticadas y cuatro recaídas) entre septiembre de 2024 y septiembre de 2025. Como esto utilizó únicamente muestras clínicas residuales desidentificadas, el comité de ética eximió el requisito de consentimiento informado por escrito.
6. Procesar inmediatamente el aspirado de médula ósea recolectado para preservar el ARN utilizando un método estándar de lisis fenol-guanidinium²¹.
   1. Transfiere la médula ósea recién recogida aspirada a un tubo de recogida que contenga un anticoagulante. Sacude el tubo para mezclar bien el aspirado y el anticoagulante.
   2. Extrae un volumen medido de la muestra anticoagulada y añádelo a un reactivo comercial de lisis fenol-guanidinio. Mantén una relación de volumen de 3 partes de reactivo de lisis por 1 pieza de muestra.
      PRECAUCIÓN: El reactivo de lisis fenol-guanidinium contiene productos químicos tóxicos y corrosivos que pueden causar quemaduras graves y daños en los tejidos. Realiza toda la manipulación de reactivos dentro de una campana extractora química mientras llevas el equipo de protección personal adecuado.
   3. Sacude vigorosamente el tubo para homogeneizar completamente la muestra y el reactivo de lisis. Asegúrate de que la mezcla esté completamente mezclada y de que no queden coágulos visibles en la solución.
   4. Congelar inmediatamente la muestra homogeneizada sumergiendo el tubo en nitrógeno líquido.
      PRECAUCIÓN: El nitrógeno líquido es extremadamente frío y puede causar una congelación severa al contacto. Lleva guantes criogénicos y una pantalla facial completa al manipular nitrógeno líquido.
7. Transfiere las muestras congeladas rápidamente a un congelador de -80 °C para su almacenamiento a largo plazo antes de la tubería posterior de aislamiento de ARN y secuenciación transcriptomática. Esto representa un punto seguro en el que el experimento puede pausarse y reiniciarse más tarde.
   NOTA: El flujo de trabajo presentado en este protocolo se centra completamente en generar firmas de resistencia computacional. No se realizó ninguna validación experimental independiente, como la PCR cuantitativa en tiempo real (RT-qPCR), sobre los genes diferencialmente expresados clave identificados a través de esta cadena específica.

2. Control de calidad de ARN y preparación de bibliotecas

1. Evalúa la integridad del ARN utilizando un sistema de electroforesis capilar microfluídico. Para este flujo de trabajo representativo, incluye muestras de ARN con un número de integridad del ARN (RIN) ≥ 6,0, una relación A260/280 entre 1,8 y 2,1, y sin picos visibles de degradación. Registra las proporciones de RIN y pureza medidas para cada muestra antes de la preparación para la biblioteca.
2. Entrada de 1 μg de ARN total por muestra para la preparación de la biblioteca. Purifica el ARNm del ARN total utilizando perlas magnéticas unidas a oligo-T para enriquecer los transcritos con cola poliA.
3. Fragmenta el ARNm enriquecido usando cationes divalentes. Incubar la mezcla a 94 °C durante 15 minutos en un buffer de reacción de síntesis de primera cadena 5X.
4. Sintetizar el ADNc de primera cadena usando cebadores hexaméricos aleatorios y una transcriptasa inversa sin actividad RNasa H.
5. Degradar la cadena plantilla de ARN usando RNasa H. Sintetizar el ADNc de segunda cadena utilizando ADN polimerasa I y dNTPs en un sistema de reacción de 20 μL.
6. Incubar la reacción de síntesis de la segunda cadena a 16 °C durante 1 hora. Centrifuga brevemente la mezcla de reacción a 2.000 x g para recoger el líquido en el fondo del tubo.
7. Convertir los salientes restantes en extremos romos utilizando actividades de exonucleasa y polimerasa. Adenilade los extremos 3' de los fragmentos de ADN y liga adaptadores con estructuras de lazos en forma de pasador para prepararse para la hibridación.
8. Purificar los fragmentos de la biblioteca usando perlas magnéticas de inmovilización reversibles en fase sólida para seleccionar preferentemente fragmentos de ADNc de 370–420 pb de longitud.
9. Realiza lavados con etanol durante la purificación de las cuentas. Centrifuga los tubos a 2.000 x g durante 30 segundos para recoger y eliminar completamente cualquier etanol residual antes de la elusión final.
10. Realiza amplificación por PCR utilizando una ADN polimerasa de alta fidelidad, cebadores universales de PCR y cebadores de índice específicos para muestra.
11. Ejecutar el perfil térmico de la PCR con una desnaturalización inicial a 98 °C durante 30 segundos. A continuación, se realizan 12 ciclos de 98 °C durante 10 s, 60 °C durante 30 s y 72 °C durante 30 s, terminando con una extensión final a 72 °C durante 5 minutos.
12. Purifica de nuevo los productos de PCR usando las perlas magnéticas. Aplica los mismos parámetros de centrifugación del paso 2.9 para obtener la biblioteca final.
13. Cuantifica la concentración inicial de la biblioteca usando un fluorómetro. Diluye la biblioteca final a una concentración de 1,5 ng/μL.
14. Mezcla bien la biblioteca diluida. Centrifugar la mezcla a 10.000 x g durante 1 minuto a 4 °C para eliminar cualquier residuo antes del análisis final.
15. Evalúa el tamaño del inserto de la biblioteca utilizando el sistema de electroforesis capilar microfluídico.
16. Cuantifique con precisión la concentración efectiva de la biblioteca mediante PCR cuantitativa en tiempo real (qRT-PCR) tras confirmar que el tamaño del inserto cumple con las expectativas. Asegúrate de que la concentración sea superior a 1,5 nM para garantizar la estabilidad de la biblioteca y la calidad de secuenciación.
    NOTA: Esto representa un punto seguro en el que el experimento puede pausarse. Las bibliotecas preparadas pueden almacenarse a -20 °C hasta el agrupamiento y secuenciación.

3. Agrupamiento y secuenciación de transcriptomas

1. Realizar el agrupamiento de las muestras codificadas por índice en un sistema automatizado de generación de clúster. Utiliza un kit comercial de cluster de extremo emparejado según las instrucciones del fabricante.
2. Secuenciar las preparaciones de la biblioteca en una plataforma de secuenciación de alto rendimiento tras la generación exitosa de clústeres. Genera 150 lecturas de pares de bases (bp) en extremos emparejados.

4. Control de calidad de datos y mapeo de lecturas

1. Evalúa la calidad de los datos en bruto (formato FASTQ) usando fastp v0.23.2 para el control de calidad y filtrado de lecturas en bruto. Registra los parámetros de la línea de comandos en un registro de análisis. En este flujo de trabajo, se generaban lecturas limpias eliminando lecturas que contenían adaptadores, lecturas con secuencias poly-N y lecturas de baja calidad usando configuraciones de filtrado idénticas entre muestras. Un comando representativo de extremo emparejado se proporciona en el Archivo Suplementario 1.
2. Procesa las lecturas en bruto a través de un software automatizado de preprocesamiento. Obtén lecturas limpias eliminando lecturas que contienen adaptadores, lecturas que contienen secuencias poly-N y lecturas de baja calidad. Utiliza parámetros de filtrado idénticos para todas las muestras y registra el número de lectura retenido, Q20, Q30 y contenido GC tras el filtrado.
3. Calcula el contenido de Q20, Q30 y GC de los datos limpios. Definir las variables potenciales del lote antes del análisis posterior, incluyendo la fecha de recogida de muestras, la fecha de extracción de ARN, el lote de preparación de la biblioteca, la línea de secuenciación y la secuenciación.
4. Evalúa los efectos por lotes mediante PCA y análisis de correlación muestra a muestra utilizando valores de expresión normalizados. Si las muestras se agrupan principalmente por variables técnicas en lugar de por estado clínico, documenta la variable afectada e inclúyela como covariable en la fórmula de diseño de expresión diferencial o aplica un método establecido de ajuste por lotes antes de la visualización posterior.
5. Adquiere el genoma de referencia (Homo sapiens, GRCh38) y los correspondientes archivos de anotación genética Ensembl release 109 para alineación de lectura.
6. Construye el índice del genoma de referencia usando HISAT2 v2.0.5.
7. Alinea las lecturas limpias de extremos emparejados con el genoma de referencia usando HISAT2 v2.0.5. Utiliza este enfoque de alineamiento consciente de la empalme para generar una base de datos de uniones de empalme basada en el archivo de anotación del modelo génico.

5. Predicción de transcripciones novedosas y cuantificación de la expresión génica

1. Ensamblar las lecturas mapeadas de cada muestra usando StringTie v1.3.3b en un enfoque basado en referencias. Utiliza esta herramienta para ensamblar y cuantificar transcritos de longitud completa que representan múltiples variantes de empalme para cada locus génico.
2. Cuenta el número de lecturas asignadas a cada gen usando featureCounts v1.5.0-p3. Utiliza la matriz de lectura entera en bruto resultante como entrada para el análisis de expresiones diferenciales posteriores.
3. Configurar featureCounts v1.5.0-p3 para datos de secuenciación de extremos emparejados usando la opción de extremos emparejados (por ejemplo, -p). Proporciona el archivo de anotación GTF descargado GRCh38 para definir los límites correctos de las características genómicas.
4. Calcula los Fragmentos por Kilobase de transcrito por Millón de lecturas mapeadas (FPKM) para cada gen. Utiliza valores FPKM solo para visualización descriptiva, PCA, visualización de mapas de calor y resúmenes exploratorios de expresiones; no se utilizan valores FPKM como matriz de entrada para las pruebas de expresión diferencial DESeq2.

6. Análisis diferencial de expresión génica

1. Realizar análisis de expresión diferencial entre los grupos recién diagnosticados y los grupos recaídos usando R v3.5.0 y el paquete DESeq2 R v1.20.0. Importar la matriz de conteo de lecturas en bruto generada en el paso 5.2 al entorno R y conservar los valores FPKM solo para visualización y análisis exploratorios.
2. Construye el objeto especializado del conjunto de datos requerido por el paquete de análisis. Ejecuta el comando específico (por ejemplo, DESeqDataSetFromMatrix()) para vincular la matriz de datos de conteo con la tabla de metadatos de muestra correspondiente.
3. Define la fórmula de diseño experimental dentro del objeto de software. Especifica el estado clínico (recién diagnosticado frente a recaído) como variable principal para la comparación (por ejemplo, diseño = ~ condición). Si una variable técnica por lotes se identifica en el paso 4.3 y no se confunde completamente con el estado clínico, inclúyela en la fórmula de diseño (por ejemplo, diseño = ~ lote + condición).
4. Ejecuta la función central de análisis diferencial de expresión (por ejemplo, DESeq()). Permitir que el software realice automáticamente la estimación del factor de tamaño, la estimación de dispersión y la prueba binomial de Wald²² negativa.
5. Extrae la tabla de resultados utilizando la función de extracción de resultados (por ejemplo, results()). Especifica el argumento del contraste para definir la comparación exacta (recaída frente a recién diagnosticado).
6. Ajusta los valores P resultantes para controlar la tasa de descubrimiento falso. Utiliza el procedimiento integrado de Benjamini y Hochberg aplicado automáticamente por el paquete de software²³.
7. Filtra la tabla de resultados extraída para aislar los genes diferencialmente expresados (DEGs) significativos. Asigna cualquier gen con un valor P ajustado < 0,05 y un logaritarítmica absoluta_{de 2} veces > 1 como significativamente expresado de forma diferencial.

7. Análisis de Enriquecimiento de la Ontología Génica (GO)

1. Realizar un análisis de enriquecimiento por Ontología Génica (GO) de los DEGs identificados usando clusterProfiler v3.8.1 y org. Hs.eg.db v3.6.0. Introduce la lista de identificadores de genes Entrez correspondientes a los DEGs significativos identificados en el paso 6.7.
2. Ejecuta la función de enriquecimiento GO (por ejemplo, enrichGO()). Especificar los parámetros requeridos, incluyendo la base de datos de organismos de fondo correspondiente (por ejemplo, OrgDb = org. Hs.eg.db), el dominio ontológico específico (Proceso Biológico, Componente Celular o Función Molecular) y el corte ajustado del valor P (0,05).
3. Asegúrese de que el algoritmo aplique las correcciones necesarias durante el cálculo del enriquecimiento. Confirma que el software corrige internamente el sesgo de longitud del gen y ajusta los valores P usando el método^{de Benjamini y Hochberg 24}.
4. Consideremos los términos GO con un valor P corregido inferior a 0,05 como significativamente enriquecidos. Genera un gráfico de puntos o un gráfico de barras utilizando las funciones de visualización integradas del paquete para mostrar los términos GO más enriquecidos.

8. Enciclopedia de Genes y Genomas de Kioto (KEGG) Análisis de Enriquecimiento de Vías

1. Utiliza un recurso de base de datos completo dedicado a comprender las funciones de alto nivel del sistema biológico para identificar las vías desreguladas. Prepara la misma lista de identificadores DEG Entrez importantes que se usan en el paso 7.1.
2. Ejecuta la función de enriquecimiento KEGG (por ejemplo, enrichKEGG()) dentro del paquete de anotación funcional R.
3. Define los parámetros críticos dentro de la llamada a la función. Establece el código del organismo estrictamente en humano (por ejemplo, organismo = 'tiene') y define el método de ajuste del valor P (por ejemplo, método pAdadjust = 'BH').
4. Extrae las vías KEGG estadísticamente significativas. Filtrar la salida para conservar solo aquellas vías que demuestren un valor P corregido menor a 0,05.
5. Visualiza las vías KEGG más enriquecidas. Utiliza las funciones integradas de trazado (por ejemplo, dotplot()) para mapear la significación estadística y el recuento de genes asociados a cada vía.

9. Análisis de Enriquecimiento de Conjuntos Génicos (GSEA)

1. Prepara la lista genética preclasificada necesaria para el análisis. Calcula la métrica de clasificación para todos los genes expresados usando el signo -log10 (valor P) multiplicado por el signo del cambio log₂ veces derivado del análisis diferencial de expresiones.
2. Lanzar una instalación local del software GSEA v4.2.3 del Broad Institute. Introduce la lista genética preclasificada recién generada en la interfaz^{de software 25}.
3. Descarga los conjuntos genéticos predefinidos necesarios. Adquiere los conjuntos de datos Gene Ontology (GO) y Kyoto Encyclopedia of Genes and Genomes (KEGG) de la Molecular Signatures Database (MSigDB, versión 7.5.1)²⁶.
4. Configura los parámetros del software para realizar la prueba estadística de enriquecimiento. Establece el número de permutaciones en 1.000 y selecciona el tipo de permutación como 'gene_set'.
5. Ejecuta el algoritmo de análisis para determinar si los conjuntos génicos predefinidos muestran una diferencia estadísticamente significativa y concordante entre los estados biológicos recién diagnosticados y los recaídos.
6. Evalúa la significación estadística de los perfiles de enriquecimiento generados. Definan conjuntos génicos significativos usando umbrales estrictos: un valor absoluto de la puntuación de enriquecimiento normalizado (NES) > 1,0, un valor P nominal < 0,05 y un valor q de tasa de descubrimiento falso (FDR) < 0,25.

10. Análisis de Redes de Interacción Proteína-Proteína (IBP)

1. Accede a la base de datos STRING para interacciones proteína-proteína conocidas y predichas. En este flujo de trabajo, el análisis PPI se realizó con STRING v11.5²⁷.
2. Introduce la lista de identificadores génicos de Entrez o símbolos génicos oficiales para los genes significativos expresados diferencialmente (identificados en el paso 6.7) en la interfaz de búsqueda de la base de datos. Selecciona a Homo sapiens como organismo objetivo.
3. Configura los parámetros de construcción de la red para asegurar que se recuperen interacciones de alta calidad. Establece la puntuación mínima requerida en la interacción a un umbral de confianza alto (puntuación > 0,700).
4. Exporta los datos de la red de interacción resultantes a un directorio local. Guarda el mapa de interacción como un archivo tabular estándar (por ejemplo, formato TSV).
5. Importa los datos de interacción exportados en Cytoscape v3.9.1 para la visualización y análisis de la^{red 28}.
6. Filtra la red construida para mejorar la claridad de la visualización y destacar los principales centros regulatorios. Eliminar cualquier nodo desconectado o gen huérfano que no muestre interacciones continuas que cumplan con el umbral de confianza establecido.

Resultados

Validación de cohortes clínicas y secuenciación

La ejecución exitosa del protocolo de extracción y preparación de bibliotecas de ARN aguas arriba (Figura 1) fue confirmada mediante métricas de rendimiento y calidad de secuenciación. En este conjunto de datos representativo, las muestras de médula ósea de cinco pacientes recién diagnosticados con LMA y cuatro pacientes con LMA recaente proporcionaron una media de aproximadamente 6,0 GB de datos en bruto por muestra. La evaluación de control de calidad (Tabla 2) confirmó que la calidad base y la profundidad de lectura cumplían con los umbrales requeridos para el análisis bioinformático posterior. Una baja integridad del ARN (por ejemplo, RIN < 6,0), bajas tasas de mapeo o un alto sesgo de degradación de transcritos representarían una calidad de entrada subóptima y comprometerían la fiabilidad del análisis diferencial de expresión posterior.

Varianza Transcriptómica Global y ACP

Para evaluar la varianza global del transcriptoma e inspeccionar la agrupación clínica, se realizó PCA sobre los datos de expresión normalizada. En este conjunto representativo, los grupos recién diagnosticados y recaídos mostraron separación en el espacio bidimensional (Figura 2A⁾²⁰, con PC1 y PC2 representando el 23,82% y el 18,75% de la varianza total, respectivamente. Los diagramas de Venn en la Figura 2B,C proporcionan un resumen descriptivo adicional de los genes detectados entre muestras dentro de los grupos recién diagnosticados y recaídos, apoyando las comprobaciones de reproducibilidad a nivel muestral antes del análisis posterior de expresión diferencial. Debido a que la cohorte era pequeña y no emparejada, la separación de PCA se interpretó como un resultado ilustrativo del flujo de trabajo en lugar de evidencia definitiva de biología específica del estado de la enfermedad.

Análisis de genes de expresión diferencial (DEG)

Aplicando los umbrales establecidos del protocolo (|_{log 2}FC| ≥ 1 y el valor P ajustado ≤ 0,05) a la salida de DESeq2 se identificaron 2.025 DEGs, que comprenden 772 genes regulados al alza y 1.253 genes a la baja en el grupo recaído (Figura 3A). Los expedientes candidatos con alta variación incluyeron FOXC1 (log₂FC = 7,55, P = 4,92 x ^10-5), HOXA11 (log₂FC = 7,76), HOXA11-AS (log₂FC = 7,23) y AXL (log₂FC = 3,50), junto con RHOB, PTX3 y CXCL8 regulados a la baja. La literatura existente vincula varios de estos genes con la relación con la respuesta de la LAM, señalización o^{terapia 13,29}; sin embargo, el flujo de trabajo actual solo los identifica como transcripciones de candidatos asociadas a recaídas. Cualquier papel mecanicista definitivo en la resistencia clínica requiere una validación funcional independiente posterior.

Enriquecimiento funcional y de vías (GO, KEGG y GSEA)

El protocolo de anotación funcional mapeaba los DEGs a sistemas biológicos más amplios. El análisis GO identificó el enriquecimiento de términos relacionados con la transducción de señales mediada por GTPasa pequeña, el transporte de iones metálicos y el ensamblaje de cromatina (Figura 4A–C). El mapeo de la vía KEGG identificó asociaciones con interacciones entre receptores ECM y citocinas-citoquinas (Figura 4D). GSEA mostró el enriquecimiento de procesos biosintéticos de ARN en el grupo recaído y el enriquecimiento de las vías del metabolismo energético en el grupo recién diagnosticado (Figura 5A). Estos resultados de enriquecimiento proporcionan una hoja de ruta descriptiva de conjuntos génicos alterados y deben interpretarse como asociaciones generadoras de hipótesis más que como factores probados de recaída.

Construcción de redes de interacción proteína-proteína (PPI)

La red STRING inicial contenía 56 nodos y 193 interacciones. Tras la eliminación de ganglios desconectados u huérfanos, la subred Cytoscape mostrada contenía 42 nodos y 136 interacciones (Figura 5B). El análisis modular de red priorizó TP53, CCL2, CXCL8 e IL6 como centros matemáticos centrales con el mayor número de interacciones. Dado que la red PPI se basa en puntuaciones de interacción predichas en bases de datos (por ejemplo, puntuación ATF3: 0,982), la identificación de los centros debe interpretarse como priorización de objetivos para futuros estudios empíricos en lugar de evidencia directa de evasión por apoptosis mediada por p53 u otros mecanismos de resistencia.

Los datos brutos de secuenciación de ARN generados en este protocolo han sido depositados en el repositorio Figshare y son accesibles públicamente a través del siguiente DOI: https://doi.org/10.6084/m9.figshare.30655814. Los datos procesados y los archivos de análisis asociados se incluyen en el artículo y sus materiales complementarios. Los parámetros representativos de línea de comandos y los ajustes de análisis usados para reproducir el flujo de trabajo computacional se proporcionan como Archivo Suplementario 1. Todos los datos que respaldan los resultados de este estudio están disponibles sin restricciones.

Identificación del paciente	Edad (años)	Sexo	Mutaciones moleculares	Supervivencia/Seguimiento (meses)	Estado clínico
R_AML_1	70	Masculino	FLT3-ITD (+)	22	Fallecido
R_AML_2	29	Femenino	NPM1 (+)	11	Viva
R_AML_3	40	Masculino	CEBPA (+)	17	Viva
R_AML_4	55	Femenino	Triple Negativo*	24	Fallecido

Tabla 1: Características demográficas y clínicas de los pacientes en el grupo de LMA recaidente (R_AML). La Tabla 1 resume las características demográficas y clínicas de la cohorte de LMA recaidente utilizada en el análisis representativo, incluyendo características clínicas a nivel de paciente relevantes para la interpretación del flujo de trabajo transcriptómico.

Ejemplo	Biblioteca	Raw_reads	Raw_bases	Clean_reads	Clean_bases	Error_rate	P20	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

Tabla 2: Resumen de la calidad de los datos. La Tabla 2 informa de métricas de calidad de secuenciación para cada muestra, incluyendo rendimiento de lectura, calidad base, contenido de GC e información de control de calidad relacionada con el mapeo utilizada para determinar si las muestras eran adecuadas para análisis posteriores.

Figura 1: Flujo de trabajo del protocolo. El flujo de trabajo resume las principales etapas experimentales y computacionales, incluyendo la recogida de muestras clínicas, control de calidad de ARN, preparación y secuenciación de bibliotecas, procesamiento y alineación de lecturas, cuantificación de transcripciones, análisis diferencial de expresiones, enriquecimiento GO/KEGG, GSEA y construcción de redes PPI. Por favor, haz clic aquí para ver una versión ampliada de esta figura.

Figura 2: Análisis cuantitativo de muestras. (A) Se realizó un análisis de componentes principales (PCA) para evaluar las diferencias entre grupos y la reproducibilidad dentro del grupo. La ACP se realizó utilizando métodos algebraicos lineales basados en valores normalizados de expresión génica en todas las muestras. (B, C) Diagramas de Venn que muestran genes detectados en muestras de los grupos AML y R_AML, respectivamente. Las regiones restringidas por muestra indican genes detectados en muestras individuales, mientras que las áreas solapadas representan genes comúnmente detectados en dos o más muestras. Por favor, haz clic aquí para ver una versión ampliada de esta figura.

Figura 3: Análisis diferencial de expresión génica. (A) Gráfico de barras que muestra el número de genes (DEGs) expresados diferencialmente entre grupos de comparación, identificados por DESeq2 con umbrales de P-valor ajustado ≤ 0,05 y |_{log 2}FoldChange| ≥ 1. (B) Trama volcánica de DEGs. El eje x representa los valores de log₂de FoldChange, y el eje y representa -log10 (valor P). Las líneas discontinuas azules indican las líneas umbral utilizadas para la selección de DEG. (C) Mapa de calor jerárquico de agrupamiento de DEGs. El eje x indica nombres de muestras, y el eje y muestra valores de expresión normalizados de los DEG. Por favor, haz clic aquí para ver una versión ampliada de esta figura.

Figura 4: Análisis de enriquecimiento funcional de genes expresados diferencialmente. (A) Parcela del bar de enriquecimiento GO. El eje x representa los términos GO, y el eje y muestra la significancia de enriquecimiento, expresada como -log10(padj). Los colores representan BP (Proceso Biológico), CC (Componente Celular) y MF (Función Molecular). (B) Parcela burbuja de enriquecimiento GO. El eje x representa la proporción de DEGs anotadas a cada término GO en relación con el número total de DEGs, y el eje y indica los términos GO. El tamaño de la burbuja corresponde al número de genes anotados, y los gradientes de color representan la significación del enriquecimiento. (C) Parcela de barra de enriquecimiento KEGG. El eje x representa las vías KEGG, y el eje y denota la significación de enriquecimiento. (D) Trama burbuja de enriquecimiento KEGG. El tamaño de la burbuja indica el número de genes anotados, y los degradados de color reflejan la importancia del enriquecimiento. Por favor, haz clic aquí para ver una versión ampliada de esta figura.

Figura 5: Enriquecimiento GSEA y análisis de redes de interacción proteína-proteína (PPI). (A) Gráfico de barras que muestra puntuaciones normalizadas de enriquecimiento (NES) para conjuntos génicos significativos seleccionados. Los valores positivos de NES indican enriquecimiento en el grupo R_AML, mientras que los valores negativos de NES indican enriquecimiento en el grupo de LMA recién diagnosticado. (B) Red de interacción proteína-proteína (IBP). Cada nodo representa una proteína, y cada arista denota una interacción entre proteínas conectadas. Por favor, haz clic aquí para ver una versión ampliada de esta figura.

Discusión

Pasos críticos en el protocolo

La ejecución exitosa de este flujo de trabajo bioinformático depende de varios pasos críticos. En primer lugar, la congelación inmediata y la correcta lisis de la aspiración de médula ósea (Paso 1.6) son fundamentales, ya que el microambiente de la médula ósea es rico en ribonucleasas que pueden degradar rápidamente la integridad^{transcriptómica 30}. Durante la fase computacional, la configuración correcta de la fórmula del diseño experimental dentro del paquete DESeq2 (Paso 6.3) es fundamental para una expresión diferencial precisa, especialmente cuando se contrasta el estado clínico (recién diagnosticado frente a recaído) mientras se consideran las posibles variables de confusión. Finalmente, aplicar umbrales rigurosos de tasa de descubrimiento de falsos (FDR) durante el Análisis de Enriquecimiento de Conjuntos Génicos (GSEA) (Paso 9.6) es un punto de control estadístico crítico para evitar la sobreinterpretación de redes funcionales falsas positivas.

Modificaciones y solución de problemas

Un desafío común en este método es la presencia de efectos por lotes, que con frecuencia ocurren cuando se recogen muestras clínicas y se secuencian a lo largo de periodos de tiempo prolongados. Las variables por lotes deben definirse antes del análisis, incluyendo la fecha de recogida de la muestra, la fecha de extracción de ARN, el lote de preparación de la biblioteca, la línea de secuenciación y la ejecución de secuenciación. Si el análisis de correlación por PCA o muestras revela agrupamiento basado en la fecha de secuenciación u otra variable técnica en lugar del fenotipo clínico, los usuarios deben modificar el protocolo incluyendo la variable por lotes en la fórmula de diseño de expresión diferencial cuando sea estadísticamente factible o aplicando algoritmos de corrección por lotes, como ComBat o SVA, antes de la visualización³¹. Si se aplica este protocolo a sangre entera en lugar de aspiraciones de médula ósea, una modificación esencial es la inclusión de un paso de agotamiento del ARNm de globina durante la preparación de la biblioteca para evitar que los transcritos de globina muy abundantes monopolicen la profundidad de lectura de secuenciación. Las versiones de software y los parámetros principales para el flujo de trabajo representativo se complementaron de la siguiente manera: 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, secuenciación de extremos emparejados a 150 pb, GSEA v4.2.3 con 1.000 permutaciones de conjuntos de genes, MSigDB v7.5.1, STRING v11.5 con interacciones de alta confianza y Cytoscape v3.9.1. Los parámetros representativos de línea de comandos y los ajustes de análisis se proporcionan en el Archivo Suplementario 1.

Limitaciones del método

Aunque es completo, este protocolo tiene limitaciones metodológicas inherentes. Primero, utiliza secuenciación de ARN a granel, que captura el perfil transcriptómico medio de toda la aspiración de médula ósea y carece de resolución espacial de célula única. Por lo tanto, el flujo de trabajo no puede determinar si una firma asociada a la recaída aumentada se origina en células madre de leucemia, células estromales, células inmunitarias o cambios en la composición del tipo celular³². Segundo, el conjunto de datos representativo es pequeño (n = 9) y no emparejado, lo que limita la robustez estadística y evita una inferencia causal definitiva. Tercero, el flujo de trabajo es puramente in silico. Genera núcleos regulatorios y vías de señalización candidatos, pero no puede validar de forma independiente su necesidad funcional en quimiorresistencia sin validación experimental ortogonal in vitro o in vivo.

Estudios genómicos recientes de células únicas y unicelulares han ampliado el marco de referencia de LMA al resolver la heterogeneidad de estado celular, la arquitectura clonal y la evolución asociada a terapias con mayor resolución 33,34,35,36. Estos enfoques son complementarios al flujo de trabajo de RNA-seq masivo descrito aquí: la secuenciación masiva proporciona una estrategia de cribado práctica y rentable para firmas transcriptómicas a nivel de cohorte, mientras que los métodos unicelulares y multiómicos pueden utilizarse en estudios de seguimiento para asignar señales candidatas a poblaciones celulares malignas o microambientales específicas.

Importancia respecto a los métodos existentes

A pesar de estas limitaciones, esta cadena transcriptómica ofrece ventajas frente a técnicas alternativas de diagnóstico y análisis. Las evaluaciones clínicas tradicionales de la recaída de LMA suelen basarse en paneles multiplex de qPCR dirigidos o citometría de flujo estándar. Aunque útiles para diagnósticos rápidos, estos métodos dirigidos están limitados por sondas predefinidas y solo pueden evaluar marcadores de resistencia^{conocidos 19}. Al utilizar secuenciación transcriptoma imparcial y genómica conjunta junto con análisis en red, este protocolo puede nominar transcritos novedosos y asociaciones a nivel de sistema que los métodos dirigidos existentes pueden pasar por alto.

Importancia y posibles aplicaciones

La metodología descrita en este protocolo es relevante para la hematología traslacional y la medicina personalizada porque puede priorizar las firmas transcriptómicas asociadas a la recaída para estudios adicionales. Una posible aplicación posterior es la nominación de antígenos superficiales o vías de evasión inmunitaria que surgen durante la recaída. Estos candidatos podrían informar el diseño de futuros estudios de validación y, si se confirman experimentalmente, podrían contribuir al desarrollo de inmunoterapias de próxima generación, incluyendo las estrategias celulares CAR-T^{o CAR-NK 37}. Estas aplicaciones traslacionales deben considerarse que generan hipótesis más que que conclusiones establecidas a partir del presente conjunto de datos.

Materiales

Lista de materiales utilizados en este artículo

Nombre	Empresa	Número de catálogo	Comentarios
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