La construcción de un modelo de nomograma pronóstico para cáncer colorrectal y la predicción de características inmunes y respuestas al tratamiento inmunológico basada en el análisis bioinformático de genes relacionados con mediadores solubles.

Según el Informe Mundial sobre el Cáncer de la Organización Mundial de la Salud, el cáncer colorrectal (CCR) se ha convertido en el tercer cáncer más prevalente y en la segunda causa principal de muertes relacionadas con el cáncer en todo el mundo, con más de 1,9 millones de nuevos casos y aproximadamente 900.000 muertes anualmente.^[1] Muchos factores pueden desencadenar el CCR, entre los cuales un estilo de vida poco saludable y los hábitos alimenticios son los principales inductores del CCR.^[2] A nivel mundial, las tasas de incidencia y mortalidad del CCR están aumentando constantemente en paralelo con el desarrollo económico y social. La Agencia Internacional para la Investigación del Cáncer predice que, para 2040, el número de nuevos casos de CCR alcanzará los 3,2 millones, y el número de muertes ascenderá a 1,6 millones.^[3] La resección quirúrgica, la terapia dirigida, la radioterapia, la quimioterapia adyuvante sistémica y la inmunoterapia son las principales opciones de tratamiento para el CCR.^[4] Actualmente, el tratamiento quirúrgico sigue siendo el método de tratamiento primario para la mayoría de los pacientes con CCR, pero los pacientes con CCR sometidos a resección quirúrgica siguen estando en riesgo de recurrencia.^[5],^[6] Por lo tanto, predecir con precisión el pronóstico del CCR para ayudar a los médicos a formular protocolos de tratamiento más personalizados se ha convertido en la clave para mejorar la supervivencia y la calidad de vida de los pacientes.

La evaluación actual del pronóstico del CCR se basa principalmente en la estadificación TNM, que es de uso común a nivel internacional; sin embargo, esta no tiene en cuenta las diferentes características patológicas de los pacientes con CCR en el mismo estadio clínico, lo que dificulta la reflexión precisa del pronóstico del paciente.^[7],^[8] Dado que el estado de los pacientes con CCR se ve influenciado por múltiples factores, la evaluación del pronóstico basada únicamente en la estadificación TNM está sesgada. Por lo tanto, se justifica una mayor innovación de un modelo de pronóstico más preciso para el CCR con el fin de mejorar la eficacia terapéutica. En los últimos años, los biomarcadores moleculares, especialmente los ARN no codificantes largos (lncARN), han demostrado un potencial significativo en el diagnóstico, la evaluación del pronóstico y el tratamiento dirigido del CCR.^[9],^[10] Ciertos lncARN (como NAMPT-AS) están significativamente sobreexpresados en los tejidos del CCR y en los exosomas del suero de los pacientes, y están estrechamente relacionados con la expresión de NAMPT. Su alta expresión se asocia con un mal pronóstico y una supervivencia reducida en el CCR.^[9] Además, múltiples lncARN participan en la metástasis hepática, la transición epitelio-mesenquimal (TEM), la respuesta inflamatoria y la resistencia a la quimioterapia/radioterapia mediante la regulación de las vías de señalización relacionadas con el tumor durante la aparición y el desarrollo del CCR. La desregulación de su expresión en los tejidos y la sangre periférica proporciona nuevas ideas para la detección no invasiva y la clasificación molecular del CCR.^[10] Estos hallazgos sugieren que los modelos de pronóstico basados en biomarcadores moleculares tienen el potencial de complementar las limitaciones de la estadificación TNM y proporcionar una estratificación de riesgos más precisa y planes de tratamiento individualizados para los pacientes con CCR.

Con el continuo avance de la investigación del cáncer, un creciente conjunto de pruebas sugiere que los mediadores solubles pueden ser factores esenciales que influyen en el pronóstico del cáncer.^[11] Una familia de pequeñas moléculas proteicas conocidas como mediadores solubles puede disolverse en los fluidos corporales y es esencial para la regulación y la mediación del sistema inmunitario. Actúan uniéndose a receptores específicos de la superficie celular, transmitiendo señales para activar o suprimir las funciones de las células inmunitarias, participando así en diversos procesos fisiológicos y patológicos, como las reacciones inflamatorias, la diferenciación celular, la proliferación y la apoptosis.^[12] Los mediadores solubles, como las citocinas, las quimiocinas y los factores de crecimiento, pueden reforzar o limitar el crecimiento de las células cancerosas al influir en la respuesta inmunitaria y la comunicación celular en el microambiente tumoral (MET).[13–15] Por ejemplo, los macrófagos asociados al tumor pueden modular la progresión de las células cancerosas mediante la liberación de mediadores solubles como CCL22, IL-10, TNF- , TGF- , IL-6 y CCL18.^[16] Las células tumorales atenúan la actividad antitumoral de las células T mediante la liberación de puntos de control inmunitarios solubles sPD-L1, facilitando así su escape inmunitario.^[17],^[18] Estudios recientes han demostrado que el pronóstico del cáncer también se ve afectado por los mediadores solubles. Por ejemplo, una alta concentración de sPD-L1 antes del tratamiento puede ser un marcador de mal pronóstico en pacientes con carcinoma de células escamosas esofágico localmente avanzado o metastásico que reciben quimioterapia citotóxica.^[19] La citocina soluble IL-8 está sobreexpresada en el CCR y está significativamente asociada con una menor supervivencia general y una menor supervivencia libre de progresión.^[20] Por lo tanto, la evaluación de los niveles de expresión de los mediadores solubles puede ser útil para la detección y la predicción del pronóstico del CCR.

Las herramientas actuales para evaluar el pronóstico del CCR siguen siendo bastante limitadas. Este estudio tiene como objetivo explorar en profundidad la importancia de los genes relacionados con los mediadores solubles (GMRS) en el pronóstico y el tratamiento del CCR, y desarrollar un modelo de evaluación del riesgo del CCR basado en los GMRS utilizando herramientas bioinformáticas, con el fin de proporcionar un nuevo método para la predicción del pronóstico. Mediante el análisis del conjunto de datos TCGA-CCR, revelamos la correlación entre el patrón de expresión de los GMRS y el pronóstico del CCR. Además, mediante la aplicación de métodos bioinformáticos como la regresión LASSO, el análisis de Cox y los algoritmos de infiltración inmunitaria, se seleccionaron 10 GMRS fiables para construir un modelo. A continuación, se verificó el modelo para su correlación con las características inmunitarias de los pacientes y los regímenes de inmunoterapia, proporcionando nuevas perspectivas para la evaluación del pronóstico del CCR.

Métodos

Fuente de datos

Este estudio utilizó los conjuntos de datos de adenocarcinoma de colon (COAD) y adenocarcinoma de recto (READ) del Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/) como conjunto de entrenamiento. Este conjunto de datos contenía datos de expresión de ARN, datos de mutación y datos clínicos de 650 muestras de tumor y 51 tejidos normales (Tabla S1). Todas las muestras tenían información completa sobre la expresión génica. En combinación con los datos de supervivencia proporcionados por TCGA, se conservaron 645 muestras que contenían información de supervivencia. Tras excluir aún más los casos con un tiempo de supervivencia ≤ 30 días, se incluyeron finalmente 599 muestras de tumor para el análisis. Los datos de recuento originales de RNA-seq se convirtieron en una matriz de expresión estandarizada de log-CPM (recuentos por millón) utilizando el paquete edgeR.

El conjunto de validación se seleccionó a partir del conjunto de datos de microarrays GSE17538 en el Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) (Tabla S2). Los criterios de selección de muestras para este conjunto de datos incluyeron: un tamaño de muestra suficiente (n = 222), información completa sobre el tiempo y el estado de supervivencia (93 muertes, 129 supervivientes), proporciones de eventos de supervivencia equilibradas y una distribución razonable del tiempo de supervivencia (0–1 año: 31 casos, 1–3 años: 60 casos, 3–5 años: 79 casos, > 5 años: 62 casos). El preprocesamiento incluyó: 1) Eliminación de los genes con una expresión media de cero; 2) Estandarización utilizando la función normalizeBetweenArrays del paquete limma. Los datos estandarizados se utilizaron como conjunto de validación externo.

Para identificar los genes candidatos relacionados con los mediadores solubles, nos referimos a métodos de investigación previos^[11] y buscamos genes relacionados con los mediadores solubles utilizando la palabra clave "mediador soluble" en la base de datos GeneCards (https://www.genecards.org.), lo que arrojó 12.018 genes característicos (Tabla S2).### Análisis diferencial de los GMRS

En el conjunto de entrenamiento, utilizamos el paquete edgeR ^[21] para realizar un análisis diferencial de los grupos normal y tumoral de CCR (la condición de selección fue |logFC| > 1, FDR < 0,05). Posteriormente, tomamos la intersección de los genes expresados diferencialmente (DEG) y los GMRS para seleccionar los GMRS significativamente diferentes, y dibujamos gráficos de volcanes y mapas de calor de expresión para mostrar sus características de expresión. A continuación, para los genes seleccionados, utilizamos el paquete clusterProfiler ^[22] para realizar análisis de enriquecimiento de Gene Ontology (GO) y Kyoto Encyclopedia of Genes and Genomes (KEGG) por separado, con el fin de explorar sus posibles funciones biológicas y vías relacionadas.### Selección de características relacionadas con el pronóstico para construir un modelo de pronóstico

Se combinaron los datos clínicos con la expresión de los DEG según las muestras, y se conservaron las muestras de tumor con un tiempo de supervivencia superior a 30 días. Se empleó el paquete survival ^[23] en el análisis de regresión univariante de los DEG. Se seleccionaron los genes candidatos en función del valor p (p < 0,001). Para evitar el sobreajuste del modelo, se utilizó el paquete glmnet ^[24] para el análisis de regresión LASSO de los genes candidatos. Empleamos la validación cruzada para seleccionar el parámetro de penalización lambda. Se eliminaron los genes con fuertes correlaciones para reducir la complejidad del modelo. El paquete survival se aplicó en el análisis de regresión multivariante de esos genes candidatos para la modelización.

Basado en este modelo, combinando los niveles de expresión de los genes candidatos y sus correspondientes coeficientes de riesgo, se calculó la puntuación de riesgo de cada paciente, y los pacientes se dividieron en un grupo de alto riesgo (HR) y un grupo de bajo riesgo (LR) en función de la mediana de la puntuación de riesgo. Las curvas de supervivencia de los dos grupos se trazaron utilizando el paquete survival. Las curvas ROC del modelo a 1 año, 3 años y 5 años se trazaron utilizando el paquete timeROC ^[25], y se calcularon los valores AUC. Finalmente, el modelo se validó externamente utilizando el conjunto de datos independiente GEO GSE17538. Se evaluó el rendimiento predictivo y la estabilidad del modelo calculando la puntuación de riesgo, trazando curvas de supervivencia y curvas ROC.### Análisis de correlación de los genes característicos

Para descubrir la relación entre los genes característicos, lanzamos un análisis de correlación de Pearson sobre los genes característicos finales. Utilizando la función rcorr del paquete Hmisc de R ^[26], calculamos los coeficientes de correlación y los valores p correspondientes entre los genes. En base a los resultados, se trazó un mapa de calor de correlación para mostrar visualmente la fuerza de las asociaciones entre los genes.### GSEA

Para determinar las funciones enriquecidas por los genes característicos, se realizó un análisis GO. En base a la puntuación de riesgo calculada por el modelo de puntuación de pronóstico, las muestras se dividieron en un grupo HR y un grupo LR. En combinación con la matriz de expresión estandarizada procesada por el paquete edgeR y la información de agrupación de riesgos, se construyó el archivo de entrada para el análisis GSEA. Posteriormente, utilizando el software GSEA, se realizó el análisis de enriquecimiento de la vía KEGG en las muestras de los dos grupos para revelar las diferencias significativas en las vías relacionadas y su posible importancia biológica.### Construcción y verificación del nomograma

Combinando las puntuaciones de riesgo y la información clínica, incluido el estadio TMN y el estadio de la enfermedad (Estadio) del CCR, trazamos un diagrama de violín para los diferentes grupos clínicos y lanzamos la prueba de Wilcoxon. Para determinar si el modelo era capaz de predecir el pronóstico de forma independiente, se combinó la información clínica y las puntuaciones de riesgo, y las muestras se sometieron a análisis de regresión univariante y multivariante para generar los correspondientes diagramas de bosque. Se aplicó el paquete rms ^[27] para trazar el nomograma para predecir las tasas de supervivencia a 1 año, 3 años y 5 años de los pacientes, y se generaron las correspondientes curvas de calibración para verificar la eficacia predictiva del nomograma.### Paisaje inmunitario de los grupos HR y LR

Para evaluar sistemáticamente las características relacionadas con la inmunidad de los pacientes con CRC, empleamos múltiples métodos para analizar los datos del transcriptoma. En primer lugar, realizamos un análisis ssGSEA en las muestras utilizando GSVA^[28] y el paquete estimate,^[29] calculamos los niveles de infiltración de células inmunitarias, la puntuación de función inmunitaria, la puntuación del estroma, la puntuación de Estimate y la pureza tumoral de los pacientes en los grupos de alto riesgo (HR) y bajo riesgo (LR), y luego trazamos los correspondientes diagramas de caja y diagramas de violín para la comparación. Se empleó el algoritmo CIBERSORT para evaluar el nivel de infiltración inmunitaria en grupos con diferentes puntuaciones de riesgo, y se representó gráficamente el diagrama de caja. Para determinar la sensibilidad de los pacientes a la inmunoterapia en diferentes grupos de riesgo, se calculó la puntuación de disfunción y exclusión tumoral (TIDE) a partir del sitio web de TIDE (http://tide.dfci.harvard.edu/), y se generaron los diagramas de violín de TIDE en los diferentes grupos. Descargamos la puntuación de inmunofenotipo (IPS) para CRC del Cancer Immunome Atlas (TCIA, https://tcia.at) y trazamos el diagrama de violín para IPS en los dos grupos.

Análisis de mutaciones somáticas

Basándonos en los datos de mutación de CRC de TCGA, se utilizó el paquete R maftools^[30] para mostrar visualmente el estado de mutación general de todas las muestras. Al mismo tiempo, se trazaron las frecuencias y los tipos de mutación de los genes candidatos en el modelo pronóstico en todas las muestras para facilitar la evaluación de las características de mutación de estos genes.

Predicción de la sensibilidad a los fármacos

A través de una serie de procesos de selección en la base de datos CellMiner (https://discover.nci.nih.gov/cellminer/), se identificaron posibles dianas farmacológicas y fármacos aplicables relacionados con la expresión de genes característicos. Específicamente, primero realizamos un análisis de correlación de Pearson para evaluar la relación entre los niveles de expresión de los genes característicos y la sensibilidad a los fármacos, y seleccionamos los fármacos con correlaciones significativas. Al mismo tiempo, basándonos en los niveles de expresión de los genes característicos, dividimos el grupo de alta expresión y el grupo de baja expresión utilizando el valor medio como umbral, y realizamos la prueba de suma de rangos de Wilcoxon para identificar los fármacos con diferencias significativas en la sensibilidad a los fármacos entre los dos grupos. Posteriormente, utilizando el paquete de software pRRophetic, utilizamos los datos de expresión génica estandarizados (log2(TPM +1)) como entrada y, basándonos en el modelo de entrenamiento de sensibilidad a los fármacos de la base de datos GDSC, predecimos los valores de IC50 de varios fármacos anticancerosos para cada muestra. Finalmente, predecimos los fármacos candidatos con valores de IC50 más bajos y una sensibilidad potencialmente mayor, proporcionando posibles referencias de tratamiento para los pacientes con CRC.

Análisis del estado de MMR en los diferentes grupos de riesgo de pacientes con CRC y su asociación con el modelo de riesgo

Para explorar aún más los beneficios de la inmunoterapia en el tratamiento de las células cancerosas en los diferentes grupos de riesgo, basándonos en el estudio previo,^[31] analizamos los niveles de expresión de los genes relacionados con MMR (MLH1, MLH2, MLH3, MLH6, PMS1, PMS2) en los pacientes con CRC de diferentes grupos de riesgo en el conjunto de datos de TCGA, las proporciones de estabilidad de microsatélites (MSS), baja inestabilidad de microsatélites (MSI-L) y alta inestabilidad de microsatélites (MSI-H), así como las diferencias en la puntuación de MSI tumoral. Al mismo tiempo, se utilizó el análisis de correlación de Pearson para evaluar la correlación entre los genes de MMR y los genes característicos del modelo pronóstico para explorar la posible relación entre la puntuación de riesgo y el estado de MMR.

Fuente de datos

Este estudio utilizó los conjuntos de datos de adenocarcinoma de colon (COAD) y adenocarcinoma rectal (READ) del Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/) como conjunto de entrenamiento. Este conjunto de datos contenía datos de expresión de ARN, datos de mutación y datos clínicos de 650 muestras de tumores y 51 tejidos normales (Tabla S1). Todas las muestras tenían información completa sobre la expresión génica. En combinación con los datos de supervivencia proporcionados por TCGA, se retuvieron 645 muestras que contenían información de supervivencia. Después de excluir aún más los casos con un tiempo de supervivencia ≤ 30 días, se incluyeron finalmente 599 muestras de tumores para el análisis. Los datos de recuento originales de RNA-seq se convirtieron en una matriz de expresión estandarizada de log-CPM (recuentos por millón) utilizando el paquete edgeR.

El conjunto de validación se seleccionó del conjunto de datos de microarrays GSE17538 en el Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) (Tabla S2). Los criterios de selección de muestras para este conjunto de datos incluyeron: tamaño de muestra suficiente (n = 222), información completa sobre el tiempo y el estado de supervivencia (93 muertes, 129 supervivientes), proporciones de eventos de supervivencia equilibradas y una distribución razonable del tiempo de supervivencia (0-1 año: 31 casos, 1-3 años: 60 casos, 3-5 años: 79 casos, > 5 años: 62 casos). El preprocesamiento incluyó: 1) Eliminación de los genes con una expresión media de cero; 2) Estandarización utilizando la función normalizeBetweenArrays del paquete limma. Los datos estandarizados se utilizaron como conjunto de validación externo.

Para identificar los genes candidatos relacionados con los mediadores solubles, nos referimos a métodos de investigación previos^[11] y buscamos genes relacionados con los mediadores solubles utilizando la palabra clave "mediador soluble" en la base de datos GeneCards (https://www.genecards.org.), lo que arrojó 12.018 genes característicos (Tabla S2).

Análisis diferencial de los SMRG

En el conjunto de entrenamiento, utilizamos el paquete edgeR^[21] para realizar un análisis diferencial de los grupos normal y tumoral de CRC (la condición de selección fue |logFC| > 1, FDR < 0,05). Posteriormente, tomamos la intersección de los genes expresados diferencialmente (DEG) y los SMRG para identificar los SMRG significativamente diferentes, y dibujamos diagramas de volcanes y mapas de calor de expresión para mostrar sus características de expresión. Luego, para los genes seleccionados, utilizamos el paquete clusterProfiler^[22] para realizar análisis de enriquecimiento de Gene Ontology (GO) y Kyoto Encyclopedia of Genes and Genomes (KEGG) por separado para explorar sus posibles funciones biológicas y vías relacionadas.

Selección de características relacionadas con el pronóstico para construir un modelo pronóstico

Se combinaron los datos clínicos con la expresión de los DEG según las muestras, y se retuvieron las muestras de tumor con un tiempo de supervivencia superior a 30 días. Se empleó el paquete survival^[23] en el análisis de regresión univariante de los DEG. Se seleccionaron los genes candidatos basándose en el valor p (p < 0,001). Para evitar el sobreajuste del modelo, se utilizó el paquete glmnet^[24] para el análisis de regresión LASSO de los genes candidatos. Empleamos la validación cruzada para seleccionar el parámetro de penalización lambda. Se eliminaron los genes con fuertes correlaciones para reducir la complejidad del modelo. El paquete survival se aplicó en el análisis de regresión multivariante de esos genes candidatos para la modelización.

Basándose en este modelo, al combinar los niveles de expresión de los genes candidatos y sus correspondientes coeficientes de riesgo, se calculó la puntuación de riesgo de cada paciente, y los pacientes se dividieron en un grupo de alto riesgo (HR) y un grupo de bajo riesgo (LR) basándose en la mediana de la puntuación de riesgo. Se trazaron las curvas de supervivencia de los dos grupos utilizando el paquete survival. Se trazaron las curvas ROC del modelo a 1 año, 3 años y 5 años utilizando el paquete timeROC^[25], y se calcularon los valores de AUC. Finalmente, el modelo se validó externamente utilizando el conjunto de datos independiente GEO GSE17538. Se evaluó el rendimiento predictivo y la estabilidad del modelo calculando la puntuación de riesgo, trazando las curvas de supervivencia y las curvas ROC.

Análisis de correlación de los genes característicos

Para descubrir la relación entre los genes característicos, lanzamos un análisis de correlación de Pearson sobre los genes característicos finales. Utilizando la función rcorr del paquete R Hmisc^[26], calculamos los coeficientes de correlación y los correspondientes valores p entre los genes. Basándonos en los resultados, se trazó un mapa de calor de correlación para mostrar visualmente la fuerza de las asociaciones entre los genes.

GSEA

Para determinar las funciones enriquecidas por los genes característicos, se realizó un análisis GO. Basándonos en la puntuación de riesgo calculada por el modelo de puntuación pronóstica, las muestras se dividieron en un grupo de HR y un grupo de LR. En combinación con la matriz de expresión estandarizada procesada por el paquete edgeR y la información de agrupación de riesgos, se construyó el archivo de entrada para el análisis GSEA. Posteriormente, utilizando el software GSEA, se realizó el análisis de enriquecimiento de la vía KEGG en las muestras de los dos grupos para revelar las diferencias significativas en las vías relacionadas y su posible significado biológico.

Construcción y verificación del nomograma

Combinando las puntuaciones de riesgo y la información clínica, incluido el estadificación TMN y el estadio de la enfermedad (Stage) del CRC, trazamos el diagrama de violín para los diferentes grupos clínicos y lanzamos la prueba de Wilcoxon. Para determinar si el modelo era capaz de predecir el pronóstico de forma independiente, se combinó la información clínica y las puntuaciones de riesgo, y las muestras se sometieron a análisis de regresión univariante y multivariante para generar los correspondientes diagramas de bosque. Se aplicó el paquete rms^[27] para trazar el nomograma para predecir las tasas de supervivencia a 1 año, 3 años y 5 años de los pacientes, y se generaron las correspondientes curvas de calibración para verificar la eficacia predictiva del nomograma.

Paisaje inmunitario de los grupos de HR y LR

Para evaluar sistemáticamente las características relacionadas con la inmunidad de los pacientes con CRC, empleamos múltiples métodos para analizar los datos del transcriptoma. En primer lugar, realizamos un análisis ssGSEA en las muestras utilizando GSVA^[28] y el paquete estimate,^[29] calculamos los niveles de infiltración de células inmunitarias, la puntuación de función inmunitaria, la puntuación del estroma, la puntuación de Estimate y la pureza tumoral de los pacientes en los grupos de HR y LR, y luego trazamos los correspondientes diagramas de caja y diagramas de violín para la comparación. Se empleó el algoritmo CIBERSORT para evaluar el nivel de infiltración inmunitaria en los grupos con diferentes puntuaciones de riesgo, y se representó gráficamente el diagrama de caja. Para determinar la sensibilidad de los pacientes a la inmunoterapia en los diferentes grupos de riesgo, se calculó la puntuación de disfunción y exclusión tumoral (TIDE) a partir del sitio web de TIDE (http://tide.dfci.harvard.edu/), y se generaron los diagramas de violín de TIDE en los diferentes grupos. Descargamos la puntuación de inmunofenotipo (IPS) para CRC del Cancer Immunome Atlas (TCIA, https://tcia.at) y trazamos el diagrama de violín para IPS en los dos grupos.

Análisis de mutaciones somáticas

Basándonos en los datos de mutación de CRC de TCGA, se utilizó el paquete R maftools^[30] para mostrar visualmente el estado de mutación general de todas las muestras. Al mismo tiempo, se trazaron las frecuencias y los tipos de mutación de los genes candidatos en el modelo pronóstico en todas las muestras para facilitar la evaluación de las características de mutación de estos genes.

Predicción de la sensibilidad a los fármacos

A través de una serie de procesos de selección en la base de datos CellMiner (https://discover.nci.nih.gov/cellminer/), se identificaron posibles dianas farmacológicas y fármacos aplicables relacionados con la expresión de genes característicos. Específicamente, primero realizamos un análisis de correlación de Pearson para evaluar la relación entre los niveles de expresión de los genes característicos y la sensibilidad a los fármacos, y seleccionamos los fármacos con correlaciones significativas. Al mismo tiempo, basándonos en los niveles de expresión de los genes característicos, dividimos el grupo de alta expresión y el grupo de baja expresión utilizando el valor medio como umbral, y realizamos la prueba de suma de rangos de Wilcoxon para identificar los fármacos con diferencias significativas en la sensibilidad a los fármacos entre los dos grupos. Posteriormente, utilizando el paquete de software pRRophetic, utilizamos los datos de expresión génica estandarizados (log2(TPM +1)) como entrada y, basándonos en el modelo de entrenamiento de sensibilidad a los fármacos de la base de datos GDSC, predecimos los valores de IC50 de varios fármacos anticancerosos para cada muestra. Finalmente, predecimos los fármacos candidatos con valores de IC50 más bajos y una sensibilidad potencialmente mayor, proporcionando posibles referencias de tratamiento para los pacientes con CRC.

Análisis del estado de MMR en diferentes grupos de riesgo de pacientes con CRC y su asociación con el modelo de riesgo

Para explorar más a fondo los beneficios de la inmunoterapia en el tratamiento de las células cancerosas en diferentes grupos de riesgo, basándonos en el estudio previo,^[31] analizamos los niveles de expresión de los genes relacionados con MMR (MLH1, MLH2, MLH3, MLH6, PMS1, PMS2) en pacientes con CRC de diferentes grupos de riesgo en el conjunto de datos TCGA, las proporciones de estabilidad de microsatélites (MSS), baja inestabilidad de microsatélites (MSI-L) y alta inestabilidad de microsatélites (MSI-H), así como las diferencias en la puntuación de MSI tumoral. Al mismo tiempo, se utilizó el análisis de correlación de Pearson para evaluar la correlación entre los genes MMR y los genes característicos del modelo pronóstico, con el fin de explorar la posible relación entre la puntuación de riesgo y el estado de MMR.

Resultados

Análisis diferencial de los SMRG

Basándonos en los datos de transcriptómica del grupo normal de CRC y el grupo tumoral en la base de datos TCGA, se obtuvieron 12.885 DEG (Tabla S4), entre los cuales 2943 DEG se superpusieron con el conjunto de SMRG (sobreexpresados: 1435; subexpresados: 1508) (Figura 1(A,B)). A continuación, investigamos las funciones biológicas y las vías de los DEG. Los resultados del enriquecimiento GO mostraron que los DEG se concentraban en gran medida en funciones como la actividad de transporte transmembrana de iones metálicos, la unión amida y la actividad endopeptidasa (Figura 1(C)). El análisis KEGG demostró que los DEG se agrupaban principalmente en la vía de interacción ligando-receptor neuroactiva (Figura 1(D)).

Construcción y verificación del modelo pronóstico

Después de combinar los datos de expresión de los DEG con los datos clínicos, realizamos un análisis de regresión univariante utilizando el paquete survival, seleccionando 16 genes candidatos relacionados con la supervivencia basándonos en un valor p < 0,001. El análisis LASSO de estos 16 genes, basado en el paquete glmnet, arrojó 15 genes característicos (Figura 2(A,B)). Al realizar un análisis de regresión multivariante sobre los 15 genes característicos, se creó un modelo basado en 10 genes característicos (Figura 2(C)). La fórmula es la siguiente:

Puntuación de riesgo = 0,147*NLGN1 + 0,1823*FJX1 + 0,4391*TIGD1 + 0,2003*ATP2A1 + 0,3226*SEZ6L2 + 0,0963*FOXD1 + 0,0729*FABP4 + 0,0932*LRP2 – 0,0589*ATOH1 + 0,0984*FTCD.

Por lo tanto, calculamos la puntuación de riesgo de las muestras clínicas, utilizamos el paquete timeROC para trazar las curvas ROC y obtuvimos los valores AUC (1 año: 0,73, 3 años: 0,77 y 5 años: 0,78) (Figura 3(A)). Se generaron dos grupos con diferentes riesgos: el grupo HR y el grupo LR, con la mediana de la puntuación de riesgo como valor crítico (Figura 3(B)). Aplicamos el paquete survival en el análisis de supervivencia y trazamos las curvas de supervivencia, revelando que, en comparación con el grupo LR, el grupo HR tenía un tiempo de supervivencia significativamente menor (Figura 3(C)). El estado de supervivencia también mostró una tendencia similar, con tasas de mortalidad más altas entre los pacientes del grupo HR en un plazo de 5 años (Figura 3(D)). Para examinar el rendimiento predictivo del modelo pronóstico, lo validamos aún más en el conjunto de datos GEO (Figura 4), revelando que los valores AUC del modelo para 1 año y 3 años fueron ambos superiores a 0,7, lo que implicó un buen rendimiento predictivo del modelo en el conjunto de validación (Figura 4(A)). El análisis de supervivencia de Kaplan-Meier (K-M) y la distribución del estado de supervivencia de los dos grupos también demostraron tendencias similares a las del conjunto de entrenamiento, con un estado de supervivencia del grupo HR peor que el del grupo LR (Figura 4(B,D)). Dados los resultados anteriores, el modelo pronóstico de CRC basado en SMRG que establecimos tiene un buen rendimiento en la predicción del pronóstico.

Análisis de correlación de los genes característicos

Un análisis adicional de los 10 genes característicos identificados en el análisis de correlación de Pearson (Figura 5) reveló que todos los pares de genes con correlaciones significativas estaban correlacionados positivamente. Entre ellos, la correlación entre el gen TIGD1 y el gen ATP2A1 fue la más fuerte (P < 0,001). El gen FJX1 se correlacionó significativamente con 7 genes. Los genes SEZ6L2 y FOXD1 se correlacionaron significativamente con 6 genes. La expresión de TIGD1 estaba estrechamente relacionada con 5 genes. La asociación de estos genes parece explicar los mecanismos biológicos, como la coexpresión, que contribuyen a los diferentes riesgos en los pacientes con CRC.

GSEA

Basándonos en el conjunto de datos TCGA-CRC, se realizó un análisis GO de los 10 genes característicos seleccionados (Figura 6(A)), revelando que estos genes estaban relacionados principalmente con funciones de unión, como la unión amida, la unión a ácidos orgánicos y la unión a ácidos carboxílicos. Los pacientes con CRC se clasificaron en grupos HR y LR según la mediana de la puntuación de riesgo, seguido de un análisis de enriquecimiento de la vía KEGG con el uso del software GSEA. Los resultados revelaron (Figura 6(B)) que el grupo HR se enriqueció significativamente en vías como la cardiomiopatía dilatada, el carcinoma de células basales, la biosíntesis de glicosaminoglicanos-condroitín sulfato y la contracción del músculo liso vascular.

Construcción y validación del nomograma

Se calcularon y analizaron las diferencias de las puntuaciones de riesgo para T, M, N y Estadio (Figura 7(A)). Los resultados demostraron que la puntuación de riesgo presentaba diferencias significativas en T, M, N y Estadio (P < 0,05), lo que indica que la puntuación de riesgo tiene el potencial de caracterizar la progresión del cáncer. Para determinar si la puntuación de riesgo podía predecir de forma independiente el pronóstico del CRC, combinamos las características clínicas (edad, T, M, N y Estadio) con la puntuación de riesgo para los análisis de regresión de Cox univariante y multivariante. Los resultados (Figura 7(B,C)) mostraron que, en ambos análisis de regresión, las razones de riesgo del estadio M y la puntuación de riesgo fueron mayores que 1 y estadísticamente significativas (P < 0,05), lo que implica que un aumento en la puntuación de riesgo condujo a un mayor riesgo de muerte en los pacientes con CRC. Esto sugirió que la puntuación de riesgo era un factor independiente importante en el CRC y prometía en la predicción del estado de supervivencia del CRC. Para proporcionar una imagen más completa del impacto de la puntuación de riesgo en la información clínica, se dibujó un nomograma del modelo para predecir la supervivencia del paciente (Figura 7(D)). Las curvas de corrección de 1 año, 3 años y 5 años demostraron que la tasa de supervivencia predicha del modelo se ajustaba bien a la tasa de supervivencia real (Figura 7(E–G)), lo que indica que el nomograma tenía un excelente rendimiento en la predicción de la supervivencia.

Análisis del paisaje inmunitario

En ssGSEA, se determinaron puntuaciones distintas de células inmunitarias y funciones inmunitarias basándose en la expresión génica en las muestras. Los resultados revelaron que, entre las 29 células inmunitarias y funciones inmunitarias analizadas, existían diferencias significativas en las puntuaciones de 9 de ellas entre el grupo HR y el grupo LR (Figura 8(A)). Las puntuaciones de las células NK y las células Th1 en el grupo LR fueron considerablemente más altas que las del grupo HR. Las puntuaciones de Macrófagos y la respuesta de tipo II de IFN en el grupo LR fueron significativamente más bajas que las del grupo HR. También se calcularon las puntuaciones de inmunidad, estroma, pureza tumoral y estimación utilizando el paquete estimate (Figura 8(B)). Los resultados mostraron que el grupo HR tenía una puntuación de estroma significativamente más alta que el grupo LR, mientras que las otras tres puntuaciones no mostraron diferencias significativas. El análisis CIBERSORT demostró que, en el grupo HR, las células T reguladoras (Tregs) y los macrófagos M0 eran considerablemente más altos que en el grupo LR, mientras que las células plasmáticas, las células dendríticas activadas y las células T de memoria en reposo eran significativamente más bajas que en el grupo LR (Figura 8(C)). El análisis de los puntos de control inmunitario mostró grandes diferencias en los niveles de expresión de 15 genes de puntos de control inmunitario entre los dos grupos. Entre ellos, CD14, ADORA2A, LAG3, LAIR1, NRP1, TNFRSF18, TNFSF4 y VTCN1 se sobreexpresaron en gran medida en el grupo HR, mientras que BTLA, CD244, CD40LG, CD44, CD48, ICOS e IDO2 se sobreexpresaron en el grupo LR (Figura 8(D)). El diagrama de caja HLA mostró niveles de expresión más altos de 2 genes HLA en el grupo LR (Figura 8(E)). Además, el diagrama de violín de las puntuaciones TIDE mostró que el grupo LR tenía puntuaciones TIDE más bajas, lo que implica una menor probabilidad de resistencia a la evasión inmunitaria tumoral, una mayor tasa de respuesta a la terapia de bloqueo de puntos de control inmunitario (ICB) y una mayor probabilidad de beneficiarse del tratamiento ICB (Figura 8(F)). El análisis estadístico de IPS se mostró mediante un diagrama de violín, que demostró que, entre los 3 IPS, el grupo LR tenía puntuaciones IPS significativamente más altas en comparación con el grupo HR, con puntuaciones IPS más altas que indican una mejor respuesta a la inmunoterapia (Figura 8(G)). En conjunto, los pacientes con CRC en el grupo LR pueden beneficiarse más de la inmunoterapia.

Análisis de la mutación somática

Se ha demostrado que las mutaciones génicas y la carga de mutación tumoral (TMB) están implicadas en la patogénesis de varios cánceres, incluido el CRC.^[32] Por lo tanto, investigamos la relación entre los genes característicos y las mutaciones génicas del CRC en el conjunto de datos TCGA utilizando el paquete maftools. En TCGA-CRC, un total de 540 muestras tumorales tenían datos de mutación. Dado que no todas las muestras tumorales en el conjunto de datos TCGA tienen datos completos de expresión de RNA-seq, información de seguimiento clínico y datos de mutación coincidentes simultáneamente, este estudio retuvo solo las muestras que tenían los tres tipos de información para garantizar la coherencia y la comparabilidad del análisis. Por lo tanto, el número de muestras efectivas para el análisis de mutación fue menor que el número total de muestras en el análisis de expresión de RNA. Al resumir las mutaciones (Figura 9(A)), encontramos que la mayoría de las mutaciones eran mutaciones sin sentido, siendo SNP el tipo de mutación más común, que se presentaba principalmente como mutaciones C>T y C>A. Un análisis de visualización adicional de las mutaciones de los 10 genes característicos en el modelo reveló que solo 128 muestras tenían resultados de detección de mutación disponibles para estos 10 genes (posiblemente debido a un secuenciamiento no cubierto, un control de calidad fallido o un secuenciamiento insuficiente). Además, el 82,03% de las muestras presentaban mutaciones génicas. Entre ellas, la frecuencia de mutación de LRP2 fue la más alta, lo que representó el 56% de las muestras mutadas (Figura 9(B)). Los diagramas de cascada de las 30 mutaciones principales demostraron que los 10 genes mutados principales eran básicamente los mismos en ambos grupos, incluidos TTN, TP53, MUC16, SYNE1, APC, FAT4, RYR2, KRAS y ? K ? (Figura 9(C,D)). Entre ellos, APC y TP53 fueron los genes mutados principales, con tasas de mutación superiores al 50%.

Análisis de la sensibilidad a los fármacos

Predijimos la relación entre los genes característicos y la sensibilidad a los fármacos basándonos en la base de datos CellMiner. Los resultados demostraron que la expresión de los fármacos Vemurafenib, PLX-4720, RAF-265 y el gen NLGN1, así como la expresión del fármaco ZM-336372 y el gen ATP2A1, estaban significativamente correlacionadas positivamente (los coeficientes de correlación fueron 0,558, 0,557, 0,431 y 0,444, respectivamente) (Figura 10(A,B)). Debido al peor pronóstico del grupo HR, utilizamos además el paquete pRRophetic para predecir la sensibilidad de tres fármacos específicos del grupo HR: PF-4708671, PI-103 y XAV939 (Figura 10(C)). La puntuación de sensibilidad estimada más alta implicó una menor sensibilidad al fármaco. Descubrimos que el grupo HR presentaba una sensibilidad significativamente mayor a los tres fármacos, lo que implicó que los tres fármacos podrían tener un potencial de aplicación en pacientes con CRC con mal pronóstico.

Análisis de las características de la inestabilidad de microsatélites (MSI), la expresión de los genes de reparación de errores de emparejamiento (MMR) y sus correlaciones en diferentes grupos de riesgo

Los estudios han demostrado que aproximadamente el 15% de los casos de CCRC presentan características de MSI, lo cual es la manifestación funcional de la deficiencia de las proteínas MMR (dMMR).^[33],^[34] La inmunoterapia suele ser eficaz en pacientes con CCRC dMMR y MSI-H.[35–37] Por lo tanto, en el conjunto de datos TCGA-CRC, comparamos los niveles de expresión de los genes MMR, la distribución de MSS/MSI-L/MSI-H y las diferencias en la puntuación de MSI entre los pacientes de diferentes grupos de riesgo. Los resultados mostraron que, entre los diversos genes MMR, solo MLH1 presentó diferencias significativas entre los grupos de alto riesgo (HR) y bajo riesgo (LR), con una mayor expresión en el grupo LR que en el grupo HR (P < 0,05, Figura 11(A)). El análisis del tipo de microsatélite indicó que ambos grupos eran principalmente MSS, con proporciones similares de MSI-H y MSI-L, y no se observó ninguna diferencia significativa (Figura 11(B)). Además, las puntuaciones de MSI de los dos grupos no mostraron diferencias estadísticas (P > 0,05, Figura 11(C)). Además, para analizar la relación entre los genes del modelo y los genes MMR, realizamos un análisis de correlación entre ellos. Los resultados revelaron que la mayoría de los genes del modelo estaban significativamente correlacionados negativamente con los genes MMR (P < 0,05, Figura 11(B)), entre los cuales los genes significativamente correlacionados negativamente con MLH1 fueron ATP2A1, FOXD1 y SEZ6L2 (P < 0,05). Además, el gen del modelo TIGD1 estaba significativamente correlacionado positivamente con los genes MMR (P < 0,05).

Análisis diferencial de los SMRG

Basándose en los datos de transcriptómica del grupo normal y del grupo tumoral de CCRC en la base de datos TCGA, se obtuvieron 12.885 DEG (Tabla S4), entre los cuales 2943 DEG se superpusieron con el conjunto de SMRG (sobreexpresados: 1435; subexpresados: 1508) (Figura 1(A,B)). A continuación, investigamos las funciones biológicas y las vías de los DEG. Los resultados del enriquecimiento GO mostraron que los DEG se concentraban en gran medida en funciones como la actividad de transporte transmembrana de iones metálicos, la unión a amidas y la actividad endopeptidasa (Figura 1(C)). El análisis KEGG demostró que los DEG se agrupaban principalmente en la vía de interacción ligando-receptor neuroactiva (Figura 1(D)).

Construcción y verificación del modelo pronóstico

Después de combinar los datos de expresión de los DEG con los datos clínicos, realizamos un análisis de regresión univariante utilizando el paquete survival, seleccionando 16 genes candidatos relacionados con la supervivencia basándonos en un valor p < 0,001. El análisis LASSO de estos 16 genes, utilizando el paquete glmnet, arrojó 15 genes característicos (Figura 2(A,B)). Al realizar un análisis de regresión multivariante sobre los 15 genes característicos, se creó un modelo basado en 10 genes característicos (Figura 2(C)). La fórmula es la siguiente:

Riskscore = 0,147*NLGN1 + 0,1823*FJX1 + 0,4391*TIGD1 + 0,2003*ATP2A1 + 0,3226*SEZ6L2 + 0,0963*FOXD1 + 0,0729*FABP4 + 0,0932*LRP2–0,0589*ATOH1 + 0,0984*FTCD.

Por lo tanto, calculamos la puntuación de riesgo de las muestras clínicas, utilizamos el paquete timeROC para trazar las curvas ROC y obtuvimos los valores AUC (1 año: 0,73, 3 años: 0,77 y 5 años: 0,78) (Figura 3(A)). Se generaron dos grupos con diferentes riesgos: el grupo HR y el grupo LR, utilizando la mediana de la puntuación de riesgo como valor crítico (Figura 3(B)). Aplicamos el paquete survival en el análisis de supervivencia y trazamos las curvas de supervivencia, revelando que, en comparación con el grupo LR, el grupo HR presentaba un tiempo de supervivencia significativamente menor (Figura 3(C)). El estado de supervivencia también mostró una tendencia similar, con tasas de mortalidad más altas entre los pacientes HR en un plazo de 5 años (Figura 3(D)). Para examinar el rendimiento predictivo del modelo pronóstico, lo validamos aún más en el conjunto de datos GEO (Figura 4), revelando que los valores AUC del modelo para 1 año y 3 años fueron ambos superiores a 0,7, lo que implicaba un buen rendimiento predictivo del modelo en el conjunto de validación (Figura 4(A)). El análisis de supervivencia de Kaplan-Meier (K-M) y la distribución del estado de supervivencia de los dos grupos también demostraron tendencias similares a las del conjunto de entrenamiento, con un estado de supervivencia peor en el grupo HR que en el grupo LR (Figura 4(B,D)). Dados los resultados anteriores, el modelo pronóstico de CCRC basado en SMRG que establecimos tiene un buen rendimiento en la predicción del pronóstico.

Análisis de correlación de los genes característicos

Un análisis adicional de los 10 genes característicos identificados en el análisis de correlación de Pearson (Figura 5) reveló que todos los pares de genes con correlaciones significativas estaban correlacionados positivamente. Entre ellos, la correlación entre el gen TIGD1 y el gen ATP2A1 fue la más fuerte (P < 0,001). El gen FJX1 estaba significativamente correlacionado con 7 genes. Los genes SEZ6L2 y FOXD1 estaban significativamente correlacionados con 6 genes. La expresión de TIGD1 estaba estrechamente relacionada con 5 genes. La asociación de estos genes parece explicar los mecanismos biológicos, como la coexpresión, que contribuyen a los diferentes riesgos en los pacientes con CCRC.

GSEA

Basándose en el conjunto de datos TCGA-CRC, se realizó un análisis GO de los 10 genes característicos seleccionados (Figura 6(A)), revelando que estos genes estaban en su mayoría relacionados con funciones de unión, como la unión a amidas, la unión a ácidos orgánicos y la unión a ácidos carboxílicos. Los pacientes con CCRC se clasificaron en grupos HR y LR según la mediana de la puntuación de riesgo, seguido de un análisis de enriquecimiento de la vía KEGG utilizando el software GSEA. Los resultados revelaron (Figura 6(B)) que el grupo HR estaba significativamente enriquecido en vías como la cardiomiopatía dilatada, el carcinoma de células basales, la biosíntesis de glicosaminoglicanos-condroitín sulfato y la contracción del músculo liso vascular.

Construcción y validación del nomograma

Se calcularon y analizaron las diferencias en las puntuaciones de riesgo para T, M, N y Estadio (Figura 7(A)). Los resultados demostraron que la puntuación de riesgo presentaba diferencias significativas en T, M, N y Estadio (P < 0,05), lo que indica que la puntuación de riesgo tiene el potencial de caracterizar la progresión del cáncer. Para determinar si la puntuación de riesgo podía predecir de forma independiente el pronóstico del CCRC, combinamos las características clínicas (edad, T, M, N y Estadio) con la puntuación de riesgo para los análisis de regresión univariante y multivariante de Cox. Los resultados (Figura 7(B,C)) mostraron que, en ambos análisis de regresión, las razones de riesgo del estadio M y la puntuación de riesgo fueron mayores que 1 y estadísticamente significativas (P < 0,05), lo que implica que un aumento en la puntuación de riesgo conduce a un mayor riesgo de muerte en los pacientes con CCRC. Esto sugiere que la puntuación de riesgo es un factor independiente importante en el CCRC y es prometedora para predecir el estado de supervivencia del CCRC. Para proporcionar una imagen más completa del impacto de la puntuación de riesgo en la información clínica, se dibujó un nomograma del modelo para predecir la supervivencia del paciente (Figura 7(D)). Las curvas de corrección de 1 año, 3 años y 5 años demostraron que la tasa de supervivencia predicha del modelo se ajustaba bien a la tasa de supervivencia real (Figura 7(E–G)), lo que indica que el nomograma tiene un excelente rendimiento en la predicción de la supervivencia.

Análisis del paisaje inmunitario

En ssGSEA, se determinaron distintas puntuaciones de células inmunitarias y funciones inmunitarias basándose en la expresión génica en las muestras. Los resultados revelaron que, entre las 29 células inmunitarias y funciones inmunitarias analizadas, existían diferencias significativas en las puntuaciones de 9 de ellas entre el grupo HR y el grupo LR (Figura 8(A)). Las puntuaciones de las células NK y las células Th1 en el grupo LR eran considerablemente más altas que en el grupo HR. Las puntuaciones de Macrófagos y la respuesta de tipo II de IFN en el grupo LR eran significativamente más bajas que en el grupo HR. También se calcularon las puntuaciones de inmunidad, estroma, pureza tumoral y estimación utilizando el paquete estimate (Figura 8(B)). Los resultados mostraron que el grupo HR tenía una puntuación de estroma significativamente más alta que el grupo LR, mientras que las otras tres puntuaciones no mostraron diferencias significativas. El análisis CIBERSORT demostró que, en el grupo HR, las células T reguladoras (Treg) y los macrófagos M0 eran considerablemente más altos que en el grupo LR, mientras que las células plasmáticas, las células dendríticas activadas y las células T de memoria CD4 eran significativamente más bajas que en el grupo LR (Figura 8(C)). El análisis de los puntos de control inmunitario mostró grandes diferencias en los niveles de expresión de 15 genes de puntos de control inmunitario entre los dos grupos. Entre ellos, CD14, ADORA2A, LAG3, LAIR1, NRP1, TNFRSF18, TNFSF4 y VTCN1 estaban muy sobreexpresados en el grupo HR, mientras que BTLA, CD244, CD40LG, CD44, CD48, ICOS e IDO2 estaban sobreexpresados en el grupo LR (Figura 8(D)). El diagrama de caja HLA mostró niveles de expresión más altos de 2 genes HLA en el grupo LR (Figura 8(E)). Además, el diagrama de violín de las puntuaciones TIDE mostró que el grupo LR tenía puntuaciones TIDE más bajas, lo que implica una menor probabilidad de resistencia a la evasión inmunitaria tumoral, una mayor tasa de respuesta a la terapia de bloqueo de puntos de control inmunitario (ICB) y una mayor probabilidad de beneficiarse del tratamiento ICB (Figura 8(F)). El análisis estadístico de IPS se mostró mediante un diagrama de violín, que demostró que, entre los 3 IPS, el grupo LR tenía puntuaciones IPS significativamente más altas en comparación con el grupo HR, con puntuaciones IPS más altas que indican una mejor respuesta a la inmunoterapia (Figura 8(G)). En conjunto, los pacientes con CCRC en el grupo LR pueden beneficiarse más de la inmunoterapia.

Análisis de la mutación somática

Se ha demostrado que las mutaciones genéticas y la carga de mutación tumoral (TMB) están implicadas en la patogénesis de diversos cánceres, incluido el CCRC.^[32] Por lo tanto, investigamos la relación entre los genes característicos y las mutaciones genéticas del CCRC en el conjunto de datos TCGA utilizando el paquete maftools. En TCGA-CRC, un total de 540 muestras tumorales presentaban datos de mutación. Dado que no todas las muestras tumorales del conjunto de datos TCGA tienen datos completos de expresión de ARN-sec, información de seguimiento clínico y datos de mutación coincidentes, este estudio retuvo solo las muestras que tenían los tres tipos de información para garantizar la coherencia y la comparabilidad del análisis. Por lo tanto, el número de muestras efectivas para el análisis de mutación fue menor que el número total de muestras en el análisis de expresión de ARN. Al resumir las mutaciones (Figura 9(A)), encontramos que la mayoría de las mutaciones eran mutaciones sin sentido, siendo la SNP el tipo de mutación más común, que se presentaba principalmente como mutaciones C>T y C>A. Un análisis de visualización adicional de las mutaciones de los 10 genes característicos en el modelo reveló que solo 128 muestras tenían resultados de detección de mutación disponibles para estos 10 genes (posiblemente debido a un secuenciamiento incompleto, un control de calidad fallido o un secuenciamiento insuficiente). Además, el 82,03% de las muestras presentaban mutaciones genéticas. Entre ellas, la frecuencia de mutación de LRP2 fue la más alta, representando el 56% de las muestras mutadas (Figura 9(B)). Los diagramas de cascada de las 30 mutaciones principales demostraron que los 10 genes mutados principales eran básicamente los mismos en ambos grupos, incluyendo TTN, TP53, MUC16, SYNE1, APC, FAT4, RYR2, KRAS y ? K ? (Figura 9(C,D)). Entre ellos, APC y TP53 fueron los genes mutados principales, con tasas de mutación superiores al 50%.

Análisis de la sensibilidad a los fármacos

Predijimos la relación entre los genes característicos y la sensibilidad a los fármacos basándonos en la base de datos CellMiner. Los resultados demostraron que la expresión de los fármacos Vemurafenib, PLX-4720, RAF-265 y el gen NLGN1, así como la expresión del fármaco ZM-336372 y el gen ATP2A1, estaban significativamente correlacionados positivamente (los coeficientes de correlación fueron 0,558, 0,557, 0,431 y 0,444, respectivamente) (Figura 10(A,B)). Debido al peor pronóstico del grupo HR, utilizamos además el paquete pRRophetic para predecir la sensibilidad de tres fármacos específicos del grupo HR: PF-4708671, PI-103 y XAV939 (Figura 10(C)). Una puntuación de sensibilidad estimada más alta implicaba una menor sensibilidad al fármaco. Descubrimos que el grupo HR presentaba una sensibilidad significativamente mayor a los tres fármacos, lo que implicaba que estos tres fármacos podrían tener un potencial de aplicación en pacientes con CCRC con mal pronóstico.

Análisis de las características de la inestabilidad de microsatélites (MSI), la expresión de los genes de reparación de errores de emparejamiento (MMR) y sus correlaciones en diferentes grupos de riesgo

Los estudios han demostrado que aproximadamente el 15% de los casos de CCRC presentan características de MSI, lo cual es la manifestación funcional de la deficiencia de proteínas MMR (dMMR).^[33],^[34] La inmunoterapia suele ser eficaz en pacientes con CCRC dMMR y MSI-H.[35–37] Por lo tanto, en el conjunto de datos TCGA-CRC, comparamos los niveles de expresión génica de MMR, la distribución de MSS/MSI-L/MSI-H y las diferencias en la puntuación de MSI entre los pacientes de diferentes grupos de riesgo. Los resultados mostraron que, entre los diversos genes MMR, solo MLH1 presentó diferencias significativas entre los grupos de alto riesgo (HR) y bajo riesgo (LR), con una mayor expresión en el grupo LR que en el grupo HR (P < 0,05, Figura 11(A)). El análisis del tipo de microsatélite indicó que ambos grupos eran principalmente MSS, con proporciones similares de MSI-H y MSI-L, y no se observó ninguna diferencia significativa (Figura 11(B)). Además, las puntuaciones de MSI de los dos grupos no mostraron diferencias estadísticas (P > 0,05, Figura 11(C)). Además, para analizar la relación entre los genes del modelo y los genes MMR, realizamos un análisis de correlación entre ellos. Los resultados revelaron que la mayoría de los genes del modelo estaban significativamente correlacionados negativamente con los genes MMR (P < 0,05, Figura 11(B)), entre los cuales los genes significativamente correlacionados negativamente con MLH1 fueron ATP2A1, FOXD1 y SEZ6L2 (P < 0,05). Además, el gen del modelo TIGD1 estaba significativamente correlacionado positivamente con los genes MMR (P < 0,05).

Discusión

Los factores de crecimiento, las quimiocinas y las citoquinas son ejemplos de mediadores solubles, que son pequeñas moléculas proteicas que controlan la inmunidad y tienen un impacto en la patogénesis de diversas enfermedades, incluido el cáncer.^[38],^[39] La concentración de mediadores solubles en la sangre y los tejidos enfermos puede influir en la activación de las vías que conducen a las respuestas inflamatorias o a la progresión de la enfermedad.^[40] Los mediadores solubles están relacionados con la regulación de la activación de las células T, la respuesta inmune humoral, la actividad de los linfocitos y el tejido estromal en el cáncer.^[17] Por lo tanto, los mediadores solubles pueden funcionar como posibles biomarcadores para diversas enfermedades, incluido el cáncer. En los últimos años, los estudios han identificado mediadores solubles como biomarcadores predictivos para el CCRC y los han aplicado en el diagnóstico precoz y el seguimiento de los pacientes con CCRC.[40–42] Sin embargo, hay pocos informes sobre modelos de predicción innovadores y precisos basados en mediadores solubles para la evaluación del riesgo de CCRC.

Basándonos en la base de datos TCGA-CRC, utilizamos los métodos LASSO y Cox para seleccionar 10 genes característicos relacionados con el pronóstico de los pacientes con CCRC a partir de 2943 genes con expresión diferencial significativa (NLGN1, TIGD1, ATP2A1, SEZ6L2, FOXD1, FABP4, LRP2, ATOH1, FTCD) y construimos un modelo de predicción del pronóstico. Se ha informado que estos genes están estrechamente relacionados con el pronóstico y la progresión tumoral de los pacientes con CCRC. Por ejemplo, FABP4 está regulado por VEGF-A y bFGF,^[43] y promueve la proliferación de las células de CCRC, el mantenimiento de las células madre y la glucólisis a través de la vía ROS/ERK/mTOR.^[44] Además, la deficiencia de expresión de FABP4 está estrechamente relacionada con la secreción de IL-6 y SerpinE1 por los adipocitos asociados al cáncer en el microambiente tumoral (TME).^[45] FOXD1 interactúa directamente con β-catenina, lo que conduce a la translocación nuclear y la activación de los genes diana posteriores (LGR5 y Sox2).^[46] Al mismo tiempo, la alta expresión de FOXD1 está significativamente correlacionada positivamente con las citoquinas relacionadas con la inmunosupresión (como IL-10, TGFB1 y su receptor TGFBR1), lo que indica que puede afectar indirectamente su expresión y función mediante la regulación de estos factores solubles. LRP2 (megaproteína) es un receptor de membrana endocítica conocido y el más grande, y su expresión anormal puede afectar el crecimiento tumoral, la metástasis y la remodelación del microambiente inmunitario mediante la regulación de la eliminación de ligandos en las células cancerosas y las células inmunitarias y estromales.^[47] En el modelo de puntuación de riesgo pronóstico relacionado con el metabolismo de los lípidos del CCRC, los pacientes con alta expresión de LRP2 tienen niveles más altos de infiltración de macrófagos asociados al tumor, fibroblastos asociados al cáncer (CAF) y células endoteliales en los tejidos tumorales, mientras que el nivel de infiltración de las células efectoras clave de la inmunidad antitumoral CD8+T es menor, lo que indica que puede promover la progresión maligna de los tumores mediante la regulación del microambiente inmunitario tumoral.^[48] La proteína de membrana de tipo I SEZ6L2 participa en el proceso de desarrollo de diversos cánceres,[49–51] y su mayor expresión en el CCRC se asocia con un mal pronóstico del paciente. La inhibición de SEZ6L2 puede desencadenar la apoptosis de las células de CCRC a través de un mecanismo dependiente de la caspasa.^[52] Además, NLGN1 promueve la penetración de las células tumorales en el monocapa endotelial mediante la regulación de la vía APC/β-catenina, mejorando así la capacidad de metástasis e invasión de las células de CCRC.^[53] Aunque los mecanismos específicos de TIGD1, ATP2A1, ATOH1 y FTCD en el CCRC son relativamente menos estudiados, los resultados del cribado de su selección como factores pronósticos independientes indican que pueden tener funciones potencialmente importantes.[54–57] Basándonos en estos hallazgos, creemos que los 10 modelos de pronóstico génico construidos en este estudio no solo tienen una capacidad predictiva estable desde el punto de vista estadístico, sino que también pueden participar en la reprogramación metabólica tumoral, la regulación del microambiente inmunitario y los procesos de proliferación y apoptosis celular a través de múltiples niveles y vías. Esto, a su vez, puede afectar conjuntamente la aparición y el desarrollo del CCRC, así como el pronóstico de los pacientes.

Mediante el análisis del modelo en el conjunto de entrenamiento (conjunto de datos TCGA) y el conjunto de validación (conjunto de datos GEO), observamos que el modelo de evaluación de riesgos demostró una alta precisión y fiabilidad en dos conjuntos de datos diferentes, lo que demuestra su potencial en las aplicaciones prácticas. Tanto los análisis de regresión de Cox multivariante como univariante demostraron además que el modelo predictivo puede funcionar como un factor independiente para evaluar eficazmente el pronóstico de los pacientes con CCRC. El nomograma, que combinó las puntuaciones de riesgo con los factores clínicos, demostró la eficacia del modelo desde una perspectiva más completa. En conclusión, estos resultados confirmaron las importantes ventajas y el valor de aplicación del modelo de evaluación de riesgos basado en SMRG en la predicción del pronóstico de los pacientes con CCRC. El TME está compuesto por componentes celulares (células estromales, células tumorales y células inmunitarias) y componentes no celulares (citoquinas y metabolitos), y es crucial para la tumorigénesis, la progresión y la resistencia a la quimioterapia.^[58],^[59] Los mediadores solubles, como reguladores de la función inmunitaria, pueden modular la progresión del cáncer al influir en el estado inmunitario del TME.[60–63] Por ejemplo, la expresión de moléculas de control inmunitario, como LAIR-1, en células inmunitarias específicas puede servir como un posible biomarcador para la progresión de los cánceres relacionados con la inflamación crónica, y también participa en la regulación del microambiente inmunitario tumoral.^[61] Además, la activación anormal de las vías de señalización y sus interacciones con los ARN no codificantes desempeñan un papel importante en la regulación inmunitaria y los mecanismos de resistencia de diversos tumores, lo que proporciona nuevos objetivos para la clasificación molecular y el tratamiento.^[62] Con el desarrollo de la medicina de precisión, la inmunoterapia individualizada y la terapia combinada basadas en las características moleculares específicas de los pacientes se han convertido en una importante dirección para mejorar la eficacia y reducir los efectos secundarios.^[63] Basándose en esto, este proyecto analizó además el paisaje inmunitario de los diferentes riesgos de CCRC utilizando ssGSEA. Descubrimos que los Tregs tenían niveles más altos de infiltración en el grupo de alto riesgo, mientras que Th1cells, T cells CD4 memory resting, NKcells, Plasma cells y Dendritic cells activated, tenían niveles más altos de infiltración en el grupo de bajo riesgo. El grupo de bajo riesgo tenía niveles más altos de infiltración inmunitaria, lo que se indica por los diferentes tipos de células entre los dos grupos. También analizamos la puntuación estromal en el TME, que es un indicador que mide los componentes no tumorales en el TME, incluida la matriz extracelular y las células inmunitarias. Descubrimos que la puntuación estromal en el grupo de bajo riesgo era considerablemente menor que en el grupo de alto riesgo, lo que apoya aún más la idea de que el grupo de bajo riesgo tenía niveles más altos de infiltración inmunitaria, lo que sugiere que los pacientes en nuestro grupo de bajo riesgo pueden tener una mejor sensibilidad a la quimioterapia y la inmunoterapia. Además, las puntuaciones TIDE e IPS para los dos grupos también confirmaron que el grupo de bajo riesgo puede beneficiarse más de la inmunoterapia.

El CCRC presenta una alta heterogeneidad a nivel genético y molecular,^[64] y esta diferencia no solo afecta a la aparición y el desarrollo de los tumores, sino que también determina la sensibilidad de los pacientes al tratamiento y su pronóstico. La tipificación molecular y la estratificación del riesgo son enfoques importantes para revelar la heterogeneidad del CCRC. Entre ellos, la prueba de MSI es uno de los indicadores de tipificación más utilizados y valiosos en la práctica clínica.^[65] El estado de MSI refleja el estado funcional del sistema MMR y tiene un valor de referencia significativo para la evaluación del pronóstico y la toma de decisiones sobre la inmunoterapia del CCRC. Los estudios han demostrado que los pacientes con dMMR/MSI-H tienen un mejor pronóstico y pueden beneficiarse significativamente de los inhibidores del punto de control inmunitario.^[37],^[65] Además, el estado de MMR no solo puede predecir el beneficio clínico de la ICB para el CCRC,^[35] sino que, en los tumores que carecen de MMR, la generación de un gran número de nuevos antígenos tumorales puede hacer que sean sensibles a este tipo de tratamiento, independientemente del origen tisular.^[36] El MSI suele estar causado por mutaciones germinales de los genes MMR (MSH2, MSH6, MLH1, PMS2) en pacientes con síndrome de Lynch o por la metilación del promotor de MLH1.^[66] Basándonos en la estrecha relación entre el MSI y el pronóstico y la respuesta al tratamiento inmunológico del CCRC, exploramos además la asociación entre los modelos de riesgo y los genes MMR y el estado de MSI. En este estudio, comparamos los niveles de expresión de los genes MMR, la distribución de MSS/MSI-L/MSI-H y las diferencias en las puntuaciones de MSI entre los diferentes grupos de riesgo. Los resultados mostraron que, entre los diversos genes MMR, solo MLH1 presentó diferencias significativas entre los grupos de alto riesgo y bajo riesgo, con una expresión significativamente menor en el grupo de alto riesgo en comparación con el grupo de bajo riesgo; sin embargo, ambos grupos eran principalmente de tipo MSS, con una proporción menor de MSI-H y MSI-L, y no hubo diferencias estadísticas en las puntuaciones de MSI. Además, el análisis de correlación reveló correlaciones significativas entre la mayoría de los genes del modelo y los genes MMR (PMS2, PMS1, MLH3, MLH1). En resumen, estos resultados sugieren que el modelo de riesgo construido en este estudio tiene una cierta asociación con la vía MMR, pero su efecto de estratificación del pronóstico no depende del estado de MSI.

Para establecer una base científica y guiar el tratamiento farmacológico del CCRC, seleccionamos algunos fármacos con conexiones significativas con los genes característicos, como Vemurafenib, PLX-4720, RAF-265 y ZM-336372. Se ha demostrado que Vemurafenib prolonga la supervivencia de los pacientes con CCRC y mejora la respuesta al tratamiento.^[67] Se ha informado que PLX4720 y ZM-336372 suprimen la mutación de BRAF.^[68] Entre ellos, el tratamiento combinado con PLX4720 y un inhibidor de PI3K puede inhibir el crecimiento de las células de CCRC con mutación de BRAF, tanto innatas como adquiridas.^[69] El fármaco RAF-265 puede reprimir la proliferación de las células de CCRC y reforzar la apoptosis mediante la supresión de la vía de señalización Raf/MEK/ERK.^[70] Considerando el rendimiento pronóstico significativamente peor de los grupos de pacientes de alto riesgo, utilizamos el paquete pRRophetic para predecir 3 fármacos sensibles al grupo de alto riesgo (PF-4708671, PI-103 y XAV939), con el fin de guiar el tratamiento farmacológico en los pacientes con CCRC.

Conclusión

En resumen, basándose en las bases de datos TCGA y GEO, este estudio analiza 10 SMRG y las utiliza para construir un nuevo modelo pronóstico para el CRC. Los pacientes con CRC pueden dividirse en grupos de alto riesgo (HR) y bajo riesgo (LR) en función del puntaje de riesgo de los 10 SMRG. La eficiencia de supervivencia del grupo de alto riesgo es significativamente peor. El posible mecanismo de los genes característicos en el CRC se revela mediante el análisis de enriquecimiento. Además, un nomograma basado en el modelo pronóstico puede predecir eficazmente las tasas de supervivencia a 1, 3 y 5 años de los pacientes con CRC. Un análisis adicional de la infiltración inmune sugirió que la inmunoterapia podría ser más beneficiosa en el grupo de bajo riesgo. Este estudio también predice posibles fármacos diana para los genes del modelo. Por lo tanto, este estudio identifica un nuevo biomarcador basado en mediadores solubles para pacientes con CRC y puede utilizarse para predecir el pronóstico.

Limitaciones

Aunque este estudio construyó un modelo de evaluación de riesgos para el CRC basado en SMRG y demostró un alto potencial en la predicción del pronóstico del paciente y la guía del tratamiento, todavía tiene varias limitaciones. En primer lugar, dado que las conclusiones del análisis se basan principalmente en muestras tumorales retrospectivas y datos de bases de datos públicas, puede haber sesgos inherentes en la selección de casos. Aunque este estudio utilizó un número suficiente de muestras de CRC como conjunto de validación para respaldar la conclusión, aún es necesario recopilar muestras clínicas a mayor escala y multicéntricas para una mayor verificación. En segundo lugar, este estudio analizó e identificó SMRG relacionados con el CRC mediante métodos bioinformáticos, pero careció de suficientes experimentos in vitro e in vivo para verificar los patrones de expresión específicos de estos genes del modelo en el CRC y su impacto directo en el pronóstico del paciente.

Perspectivas futuras

La investigación futura debe centrarse en varias áreas clave: (1) recopilar muestras clínicas prospectivas y multicéntricas a mayor escala para validar la aplicabilidad y la estabilidad de nuestro modelo en diversas poblaciones; (2) realizar experimentos in vitro e in vivo para dilucidar el papel mecanístico de los SMRG en la progresión del CRC y la regulación del microambiente inmunitario; y (3) integrar conjuntos de datos multiómicos, incluidos transcriptómica, proteómica y metabolómica, para optimizar la precisión predictiva y explorar el potencial de traslación del modelo para la terapia personalizada.

Además, desde una perspectiva de sostenibilidad, la reutilización de fármacos ofrece una estrategia prometedora para mejorar el manejo del CRC. Por ejemplo, las terapias basadas en GLP-1 o los inhibidores de la DPP-4 han demostrado que modulan la homeostasis del colesterol a través de la regulación mediada por LXR de la expresión de ABCA1/ABCG1 en los adipocitos, mejorando así los perfiles metabólicos con posibles beneficios anticancerígenos. ^[71] Dirigirse al proteasoma, en particular al núcleo catalítico 20S, utilizando inhibidores del proteasoma aprobados clínicamente o nuevos, representa otra vía terapéutica para el control tumoral. ^[72] Además, los compuestos naturales como la prodigiosina, que exhiben efectos moduladores del TME y inmunorreguladores con baja toxicidad para las células normales, tienen un potencial considerable como terapias adyuvantes. ^[15] De manera similar, el hinokitiol ha demostrado la capacidad de inhibir la progresión de la capacidad de formación de tumores, inducir la apoptosis e impedir el flujo autofágico en las células de cáncer de mama, lo que sugiere su valor profiláctico y terapéutico. ^[73] La incorporación de estos enfoques podría ampliar el alcance de la terapia del CRC al tiempo que mejora su sostenibilidad e impacto traslacional.