$$\rightleftharpoonup{xx}$$
$$\longleftharp{xx}$$,
$$\longrightharp{xx}$$,
Recuperación y preprocesamiento de datos
Este estudio se realizó utilizando datos transcriptómicos y genéticos disponibles públicamente; No se implicó directamente en humanos ni animales. Se recuperaron conjuntos de datos transcriptómicos relevantes para la resistencia al tratamiento del cáncer de mama HER2+ de la base de datos NCBI Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/)11. Se seleccionaron dos conjuntos de datos de RNA-seq, GSE231524 y GSE231525, debido a su enfoque específico en la resistencia impulsada por HER3 y la inhibición de DUSP6 en líneas celulares de cáncer de mama HER2+ (BT474 y MDA-MB-453). Estos conjuntos de datos incluían fenotipos parentales, tolerantes y resistentes a fármacos derivados de la exposición a Lapatinib (1 μM) y el knockdown a DUSP6. Las matrices de recuento en bruto y los archivos de metadatos correspondientes se accedían usando los paquetes GEOquery (v2.70.0) y Biobase (v2.62.0) en RStudio (v4.3.2)12. Los metadatos se seleccionaron para definir dos contrastes principales para cada conjunto de datos: GSE231524 compararon el control (BT474 parental, Día 0) con muestras tolerantes y resistentes a fármacos (Día 9–Mes 9), mientras que GSE231525 compararon el control (Scrambled siRNA) con DUSP6 knockdown (DUSP6-KD). El control de calidad y la normalización de datos se realizaron utilizando el marco DESeq2 (v1.42.0), que aplica una transformación estabilizadora de varianzas (VST) para reducir la heteroscedasticidad y garantizar la comparabilidad entre muestras. La distribución de datos y los patrones de agrupamiento se evaluaron visualmente utilizando ggplot2 (v3.5.0) y pheatmap (v1.0.12) para confirmar la uniformidad de los datos e identificar posibles valores atípicos antes del análisis diferencialde expresión 13,14.
En este estudio, se mantuvo una clara distinción entre los hallazgos derivados de conjuntos de datos transcriptómicos de líneas celulares y los obtenidos de conjuntos clínicos derivados de pacientes. Los datos de líneas celulares se utilizaron principalmente para análisis exploratorios, incluyendo la identificación de genes expresados diferencialmente y la generación de conocimientos mecanicistas preliminares en modelos experimentales controlados. En cambio, se emplearon conjuntos de datos derivados de los pacientes para la validación externa de patrones de expresión génica y la evaluación de la relevancia clínica, incluyendo la evaluación pronóstica. En consecuencia, los resultados de los modelos de líneas celulares y cohortes clínicas se interpretan por separado para evitar generalizaciones excesivas y garantizar un contexto traduccional adecuado para todos los hallazgos.
Análisis diferencial de expresión génica
Se realizó un análisis de expresión diferencial para identificar genes que estaban significativamente modulados entre las condiciones de control y de tratamiento. Los recuentos normalizados se procesaron utilizando el modelo de regresión ajustada (ARM) integrado en DESeq2 para estimar con precisión los cambios de pliegue log₂ y la significación estadística. La fórmula de diseño se definió como ~condición, que representa los grupos de control frente a los grupos tratados. Se consideraron que los genes con un valor p ajustado (FDR) < 0,05 y un cambio absoluto de fold log₂ ≥ 1 se expresaron de manera significativamente diferencial. La contracción por cambio de pliegue log₂ se realizó utilizando el método apeglm para mejorar la robustez en la estimación del tamaño del efecto. Los resultados del análisis se visualizaron utilizando EnhancedVolcano (v1.22.0)15 y ggplot216, que generaron gráficos de volcanes y gráficos de MA mostrando la relación entre la magnitud de expresión y la confianza estadística. También se evaluaron estimaciones de dispersión dentro de DESeq2 para asegurar un modelado preciso de varianzas y una normalización consistente entre réplicasbiológicas 17.
Recuperación e identificación de genes mitocondriales relacionados con el estrés oxidativo y expresados diferencialmente (MOS-DEGs)
Para investigar la relación entre el metabolismo energético, el estrés oxidativo y la resistencia a fármacos, se compiló una lista exhaustiva de genes mitocondriales y asociados al estrés oxidativo a partir de múltiples bases de datos, incluyendo Human MitoCarta3.018 (https://personal.broadinstitute.org/scalvo/MitoCarta3.0/human.mitocarta3.0.html), Gene Ontology (GO:0006979, response to oxidative stress) (http://geneontology.org/), la vía de fosforilación oxidativa de la Kyoto Encyclopedia of Genes and Genomes (KEGG) (https://www.genome.jp/kegg/), y la Base de Datos de Genes de Estrés Oxidativo Humano (HOSGDB) (http://hosgdb.com/). Todos los genes recuperados se estandarizaron con símbolos génicos aprobados por HGNC usando org. Hs.eg.db (v3.18.0) y AnnotationDbi (v1.64.0), mientras que se eliminaron entradas duplicadas, pseudogenes y ARN no codificantes para garantizar la precisión de las anotaciones. El panel de genes de estrés oxidativo mitocondrial (genes MOS) resultante curado se utilizó posteriormente como conjunto de referencia para la integración con los genes diferencialmente expresados identificados en ambos conjuntos de datos transcriptómicos.
La intersección de la lista de genes MOS seleccionada con los DEGs obtenidos de GSE231524 y GSE231525 se realizó en R usando funciones intersect() dplyr (v1.1.3)19 y base R. Este enfoque integrador permitió la identificación de MOS-DEGs, que representan genes funcionalmente vinculados al metabolismo mitocondrial, la regulación redox y la adaptación al estrés oxidativo. La superposición entre conjuntos de datos se visualizó utilizando el paquete VennDiagram (v1.7.3) en R para ilustrar genes compartidos y únicos entre los modelosexperimentales 20. La lista refinada de MOS-DEGs se utilizó para análisis posteriores, proporcionando una visión mecanicista sobre la reprogramación transcripcional y metabólica subyacente a la resistencia terapéutica dirigida a HER2.
Perfilado de expresiones y visualización de MOS-DEGs
El perfilado de expresión de los MOS-DEGs identificados se realizó utilizando los encapsulados ComplexHeatmap (v2.18.0)21 y pheatmap (v1.0.12) en RStudio para visualizar patrones de expresión globales en condiciones parentales, tolerantes a fármacos y resistentes. Los datos normalizados de recuento se transformaron mediante escalado de puntuación z para estandarizar la matriz de expresión génica entre muestras. El agrupamiento se realizó utilizando métodos de distancia euclidiana y de enlace completo para detectar patrones de coexpresión y distinguir perfiles transcripcionales específicos de cada condición. Se generaron mapas de calor y diagramas de agrupamiento usando ggplot2 para asegurar una diferenciación visual clara entre condiciones. Este enfoque de visualización facilitó la identificación de grupos génicos asociados a la actividad mitocondrial, la modulación del estrés oxidativo y la reprogramación metabólica bajo estados resistentes a fármacos.
Enriquecimiento funcional y anotación de vías
Para explorar la importancia biológica y los mecanismos regulatorios de los MOS-DEG identificados, se realizaron análisis de enriquecimiento de la Ontología Genética (GO) y la Enciclopedia de Genes y Genomas de Kioto (KEGG) utilizando R Studio (versión 4.3.1). Los análisis se realizaron en el entorno tidyverse utilizando múltiples paquetes Bioconductor para computación y visualización reproducibles. La anotación genética y el mapeo de identificadores se realizaban usando la organización. Hs.eg.db base de datos (https://bioconductor.org/packages/org.Hs.eg.db/) basada en el genoma de referencia Homo sapiens (GRCh38). El análisis de enriquecimiento GO se realizó utilizando el paquete clusterProfiler (versión 4.8.1; https://bioconductor.org/packages/clusterProfiler/), que categoriza los genes en tres ontologías principales: Proceso Biológico (BP), Componente Celular (CC) y Función Molecular (MF). La funciónenrichGO 22 se utilizó con parámetros configurados a valor p < 0,05 y valor p ajustado. (FDR) < 0,05, aplicando el método de corrección de Benjamini–Hochberg. Las visualizaciones, incluyendo diagramas de barras, diagramas de puntos y diagramas de cuerdas, se generaron utilizando enrichplot (https://bioconductor.org/packages/enrichplot/), ggplot213 (https://cran.r-project.org/web/packages/ggplot2/) y GOplot (https://cran.r-project.org/web/packages/GOplot/). Estas herramientas proporcionaron una visión estructurada de los términos enriquecidos de GO y sus asociaciones génicas.
El enriquecimiento de vías KEGG se realizó utilizando la función enrichKEGG() dentro del paquete clusterProfer, referenciando la base de datos humana KEGG (https://www.genome.jp/kegg/). El paquete KEGGREST (https://bioconductor.org/packages/KEGGREST/) se utilizaba para la recuperación y anotación de datos de caminos. Se consideraron significativas vías con un valor p ajustado. (valor q) < 0,05. La visualización y el mapeo de rutas se realizaron usando pathview (https://bioconductor.org/packages/pathview/), ggplot2 y enrichplot, mientras que igraph y ggraph se usaron para la representaciónde la red 23. Todos los análisis y visualizaciones de enriquecimiento se implementaron en R Studio (v4.3.1) utilizando código reproducible y flujos de trabajo estandarizados de Bioconductor, asegurando la identificación fiable de categorías funcionales enriquecidas y vías biológicas asociadas a los MOS-DEGs.
Validación basada en ROC de biomarcadores predictivos en cáncer de mama
Para validar el poder predictivo clínico de los MOS-DEGs, se realizó un análisis de curvas de características de funcionamiento del receptor (ROC) utilizando la herramienta online ROCplotter (https://www.rocplot.org/)24. ROCplotter es una plataforma web integrada que combina datos de expresión génica con conjuntos de datos clínicamente anotados de respuesta al tratamiento de 3.104 pacientes con cáncer de mama, incluyendo aquellas tratadas con quimioterapia, terapia hormonal o agentes anti-HER2.
El análisis se realizó utilizando los parámetros "respuesta completa patológica" como variable de resultado y "cualquier quimioterapia" como categoría de tratamiento. Los valores de expresión génica derivados de los conjuntos de datos de microarrays de Affymetrix se estratificaron automáticamente en grupos de respondedores y no respondedores basándose en anotaciones clínicas dentro de la plataforma.
Se aplicaron la curva de característica operativa del receptor (ROC) (AUC), la prueba U de Mann–Whitney, el cambio de pliegue y la prueba de chi-cuadrado para evaluar la capacidad de cada gen para distinguir a los respondedores de los no respondedores. El área bajo la curva (AUC) se utilizó como métrica principal para evaluar el rendimiento discriminativo. Se consideraron significativos valores de AUC superiores a 0,55 con valores p ROC < 0,05, representando un rendimiento predictivo moderado típico de biomarcadores transcriptómicos, mientras que se aplicó corrección de tasa de descubrimiento de falsos (FDR) para mantener el rigor analítico.
Todos los MOS-DEG seleccionados fueron consultados usando sus correspondientes identificadores de sonda Affymetrix. El potencial discriminativo de cada gen se evaluó entre cohortes clínicas de cáncer de mama, donde los datos de expresión se estratificaron en grupos de respondedores y no respondedores. Las curvas ROC, diagramas de caja y resultados estadísticos asociados se generaban directamente por la plataforma ROCplotter y se exportaban para visualización y comparación posteriores. El análisis cuantificó el valor predictivo de los reguladores redox y metabólicos implicados en el estrés oxidativo mitocondrial. Los genes que mostraron significación predictiva consistente entre muestras clínicas se conservaron para su inclusión en el panel predictivo final.
Análisis de expresión diferencial en tejidos tumorales, normales y metastásicos (análisis TNMplot)
Los patrones de expresión de los MOS-DEG superiores se analizaron en tejidos mamarios normales, tumorales y metastásicos utilizando la herramienta web TNMplot v2 (https://tnmplot.com/analysis/)25. Se examinaron tanto conjuntos de datos de RNA-Seq (TCGA + GTEx + MET500) como de chips genéticos para garantizar la validación multiplataforma. El módulo "Análisis Múltiple de Genes" se utilizó junto con el carcinoma invasivo de mama como el tipode tejido seleccionado 26. Los valores de expresión se transformaron log₂ y se compararon entre los grupos Tumor vs. Normal (TvsN), Metastásico vs. Tumor (MvsT) y Metastásico vs. Normal (MvsN). El TNMplot calculó automáticamente el cambio de plegamiento (FC) y los valores p utilizando la prueba U de Mann–Whitney para evaluar la significación estadística. Las distribuciones de expresión se visualizaron como diagramas de caja y diagramas de densidad generados directamente desde la interfaz del diagrama TNM, con verde, rojo y gris representando tejidos normales, tumorales y metastásicos, respectivamente. Todas las figuras se exportaron en alta resolución para su integración en la sección de resultados. Este análisis de doble plataforma permitió la identificación y validación robusta de reguladores metabólicos mitocondriales, que se asocian con la progresión del cáncerde mama 27.
Supervivencia y análisis pronóstico usando el trazador de Kaplan–Meier
Para evaluar la relevancia pronóstica de los MOS-DEG en el cáncer de mama, se realizó un análisis de supervivencia utilizando la herramienta online Kaplan–Meier Plotter (https://kmplot.com/analysis/)28. Esta base de datos integra datos de expresión génica y supervivencia de más de 4.900 pacientes con cáncer de mama, derivados de múltiples conjuntos de datos GEO, EGA y TCGA. El análisis se realizó para la supervivencia sin recurrencia (RFS) utilizando identificadores individuales de la sonda Affymetrix correspondientes a los genes priorizados: 225609_at (GSR), 201761_at (MTHFD2), 201619_at (PRDX3/AOP1) y 201128_s_at (ACLY). Los pacientes se dividieron en grupos de alta y baja expresión según el corte mediano de expresión, y se estimaron las probabilidades de supervivencia utilizando el método de Kaplan–Meier. La prueba logarítmica se utilizó para evaluar la significación estadística entre curvas de supervivencia, y la herramienta calculó automáticamente razones de riesgo (HR) con intervalos de confianza (IC) del 95%. Todos los análisis se realizaron utilizando el endpoint del RFS, sin restricción basada en el receptor hormonal o el estado HER2 (ER, PR, HER2 = todos). Se eliminaron muestras redundantes y se verificaron las suposiciones de riesgos proporcionales para garantizar la robustez estadística. Los filtros de control de calidad excluyeron microarrays sesgados. No se aplicó selección manual de la sonda ni corrección de valores p para múltiples pruebas, de acuerdo con los ajustes predeterminados del KM Plotter. La significación estadística se definió como p. < 0,05. Se visualizaron y descargaron gráficos de supervivencia en alta resolución para una interpretación más detallada, comparando los resultados entre los expresores altos y bajos de cada gen candidatoMOS 29,30.
El presente análisis se inició utilizando conjuntos de datos que comprenden exclusivamente muestras de cáncer de mama HER2+ para la identificación de genes diferencialmente expresados (DEGs) y genes hub. Posteriormente, se realizó un análisis de supervivencia sin restricciones al estado HER2 (ER, PR, HER2 = todos) para evaluar la relevancia pronóstica más amplia y la generalizabilidad de los genes identificados. Este enfoque se empleó como un paso secundario de validación en lugar de redefinir el enfoque del estudio. Por lo tanto, las implicaciones pronósticas de los genes hub identificados se interpretan con cautela, manteniendo las conclusiones primarias específicas para el cáncer de mama HER2+ .
Las secuencias de transcripción canónicas de MTHFD2-201 (ENST00000394053.7) y PRDX3-201 (ENST00000298510.4) se recuperaron del Ensembl Genome Browser (https://www.ensembl.org)31,32. La anotación y clasificación de variantes se realizó utilizando el Predictor de Efecto Variante Ensembl (VEP) (https://www.ensembl.org/vep), que proporcionó un contexto genómico detallado, alteraciones de codones y sustituciones de aminoácidos para cada variante identificada. Solo se seleccionaron variantes de error de sentido (SNPs no sinónimos) para el análisis posterior.
Predicción de la patogenicidad y priorización de variantes
Las consecuencias funcionales de cada nsSNP se evaluaron utilizando una combinación de herramientas de predicción computacional. Se aplicó SIFT (https://sift.bii.a-star.edu.sg) para evaluar la conservación de aminoácidos, clasificando variantes con una puntuación ≤ 0,05 comoperjudiciales 33. PolyPhen-2 (http://genetics.bwh.harvard.edu/pph2) estimó el impacto estructural y evolutivo de las sustituciones, donde puntuaciones ≥ 0,85 indicaban dañoprobable 34. CADD (https://cadd.gs.washington.edu) proporcionó una puntuación compuesta de pernocidad que integra múltiples anotaciones, con valores ≥ 20 que denotan un alto potencial patógeno35. Se integraron métricas complementarias de MetaLR36, Mutation Assessor y REVEL desde la interfaz VEP para mejorar la fiabilidad de lapredicción 37. Las variantes que cumplían los umbrales MetaLR ≥ 0,70, Mutation Evaluador ≥ 3,5 y REVEL ≥ 0,75 se priorizaron como probablemente patógenas.
Predicción del impacto estructural y mecanicista
Para evaluar cómo las sustituciones de aminoácidos influyen en la integridad estructural y la función bioquímica, cada nsSNP mejor clasificado fue analizado a fondo utilizando MutPred2 (http://mutpred.mutdb.org)38 y DynaMut (http://biosig.unimelb.edu.au/dynamut)39. MutPred2 estimó la probabilidad de alteración funcional, incluyendo alteración de la actividad catalítica, ganancia o pérdida de residuos de unión a metales, cambios en la accesibilidad del disolvente y modulación alostérica, con puntuaciones ≥ 0,80 clasificadas como altamente patógenas. DynaMut calculó el cambio de energía libre de Gibbs (ΔΔG) entre proteínas de tipo salvaje y mutantes, evaluando la dirección y magnitud de la alteración de la estabilidad, y generó visualizaciones de desplazamientos atómicos y reordenamientos de enlaces de hidrógeno.
Modelado estructural secundario y 3D y perfilado de accesibilidad con disolventes
Las estructuras cristalinas resueltas experimentalmente de MTHFD2 y PRDX3 se obtuvieron del Banco de Datos de Proteínas (PDB) y se procesaron utilizando PyMOL V:3.1 (https://pymol.org)40 para visualizar la distribución espacial de residuos perjudiciales. Los modelos mutantes se crearon introduciendo las sustituciones correspondientes de aminoácidos, seguidas de refinamiento estructural y minimización de energía. La inspección comparativa 3D puso de manifiesto cambios en los elementos secundarios, contactos interatómicos alterados y la proximidad espacial de los nsSNP a dominios catalíticos y de unión a cofactores, revelando posibles alteraciones en las funciones redox y metabólicas.
Se realizaron análisis secundarios de estructura y exposición a disolventes utilizando PSIPRED V: 3.2 (http://bioinf.cs.ucl.ac.uk/psipred)41,42 y NetSurfP 3.0 (https://services.healthtech.dtu.dk/service.php?NetSurfP-2.0)43. Estas herramientas predecían α-hélices, hebras β, bobinas y regiones desordenadas junto con puntuaciones relativas de accesibilidad al disolvente (RSA). Se mapearon los residuos que presentaban valores de RSA moderados a altos y un orden estructural para identificar posiciones expuestas al disolvente y funcionalmente críticas. Los sitios afectados se visualizaron en diagramas topológicos 2D para determinar si ocurrían mutaciones perjudiciales en núcleos catalíticos rígidos o en regiones de bucles flexibles, prediciendo así sus probables efectos en la dinámica del plegamiento de proteínas y la eficiencia enzimática.
Validación cruzada de bases de datos, integración funcional y validación de estabilidad
Cada nsSNP priorizado fue cruzado con bases de datos genómicas a nivel poblacional, incluyendo dbSNP, 1000 Genomes, ExAC y gnomAD, para confirmar la frecuencia de variantes, la distribución global de alélicos y asociaciones clínicas previamente reportadas. La integración de la conservación evolutiva, el modelado estructural y la predicción funcional basada en aprendizaje automático permitió identificar variantes perjudiciales de alta confianza en MTHFD2 y PRDX3. Estas mutaciones de alto impacto fueron posteriormente mapeadas en dominios funcionales para elucidar su posible papel en el desequilibrio del estrés oxidativo mitocondrial, la alteración de la señalización metabólica y la resistencia terapéutica en el cáncer de mama. Para confirmar aún más las consecuencias termodinámicas de cada sustitución perjudicial, se utilizó iMutant 3.0 (https://folding.biofold.org/i-mutant/i-mutant3.0.html)44 para predecir los efectos de las mutaciones en la estabilidad de las proteínas mediante datos de secuencia y estructuras. El análisis calculó valores de ΔΔG (kcal/mol) que representan el cambio en la energía libre entre proteínas de tipo salvaje y mutantes. Las variantes que mostraban valores negativos de ΔΔG se clasificaron como mutaciones desestabilizadoras, indicando una reducción de la estabilidad de la proteína y un aumento de la probabilidad de despliegue. La integración de las predicciones de iMutant con los resultados de DynaMut y MutPred2 proporcionó validación cruzada para identificar residuos estructuralmente críticos que probablemente afecten a la función redox, la integridad catalítica y la estabilidad conformacional general de proteínas.
Interpretación funcional integrada y relevancia terapéutica
Todos los NSSNP perjudiciales identificados fueron validados mediante referencias cruzadas con las bases de datos poblacionales dbSNP, gnomAD y ExAC para verificar frecuencias alelares menores y asociaciones previamente reportadas con fenotipos cancerosos. La interpretación integrativa de datos de conservación evolutiva, modelado estructural y estabilidad indicó que las mutaciones de alto impacto rs1471336772 (MTHFD2) y rs747786383 (PRDX3) ejercen los efectos perjudiciales más fuertes sobre la conformación proteica y la eficiencia catalítica. Los resultados computacionales sugieren colectivamente que las mutaciones en MTHFD2 desestabilizan el metabolismo redox dependiente de NADPH, mientras que las mutaciones en PRDX3 afectan la defensa contra el estrés oxidativo mediada por la peroxidasa, contribuyendo a la disfunción mitocondrial y a la agresividad tumoral. Este análisis estructural y funcional basado en nsSNP proporciona una base computacional para futuras pruebas terapéuticas y validación mutacional, destacando MTHFD2 y PRDX3 como biomarcadores de precisión para la terapia contra el cáncer de mama dirigida a redox. Para mejorar la claridad y ofrecer una visión completa de la estrategia analítica, en la Figura 2 se presenta un flujo de trabajo esquemático que resume los principales pasos del estudio. El flujo de trabajo integra análisis diferencial de expresión génica, filtrado génico mitocondrial, construcción de redes de interacción proteína-proteína, validación clínica mediante análisis ROC y caracterización estructural basada en nsSNP. Este marco paso a paso destaca la progresión lógica desde el procesamiento de datos transcriptómicos hasta la identificación de biomarcadores e interpretación funcional.
Figura 2. Flujo de trabajo integrativo de varios pasos para la identificación y validación de biomarcadores relacionados con el estrés oxidativo mitocondrial en el cáncer de mama HER2+ . Este esquema resume la línea analítica empleada en el estudio. En primer lugar, se realizó un análisis diferencial de expresión génica (DEG) en conjuntos de datos de RNA-seq (GSE231524 y GSE231525) para identificar genes significativamente alterados. Estos DEGs se intersectaron con genes relacionados con el estrés oxidativo mitocondrial curados para obtener MOS-DEGs. A continuación, se realizó un análisis de redes de interacción proteína-proteína (PPI) utilizando STRING y Cytoscape para identificar genes central y módulos funcionales. Posteriormente, se aplicó el análisis de la curva de características de funcionamiento del receptor (ROC) utilizando la plataforma ROCplotter para evaluar el rendimiento predictivo de genes seleccionados en cohortes clínicas. Finalmente, se realizaron análisis de SNP no sinónimos (nsSNP) y modelización estructural para evaluar el posible impacto funcional y estructural de variantes clave en genes priorizados (MTHFD2 y PRDX3). Este flujo de trabajo integrativo vincula análisis transcriptómicos, de red, clínicos y estructurales para identificar posibles biomarcadores y objetivos terapéuticos. Por favor, haz clic aquí para ver una versión ampliada de esta figura.
