Declaración del comité de revisión institucional
Este estudio se realizó de acuerdo con la Declaración de Helsinki. El protocolo fue aprobado por el Comité de Ética del Hospital Luohu de Medicina Tradicional China de Shenzhen (número de aprobación: 2024-LHQZYYYXLL-KY-039), y se obtuvo el consentimiento informado por escrito de todos los participantes antes de su inclusión. Los detalles de las herramientas y materiales de investigación utilizados en este protocolo se proporcionan en la Tabla de Materiales.
Origen y procesamiento de los datos
Los conjuntos de datos de expresión génica relacionados con la EPOC se obtuvieron del Gene Expression Omnibus (GEO). El conjunto de datos GSE54837 se utilizó como el conjunto de datos del transcriptoma, y el conjunto de datos GSE112811 sirvió como conjunto de validación (Tabla 1). Los genes relacionados con ac4C (ac4C-RGs) se recopilaron de la literatura18. Los genes diferencialmente expresados (DEG) entre los grupos de EPOC y control se identificaron utilizando el paquete limma de R. Se consideró que los DEG eran estadísticamente significativos si |log2FC| > 0 y p < 0,05. Se generaron gráficos de volcan para visualizar la distribución general de los cambios en la expresión génica.
Construcción de WGCNA
Se realizó un análisis de WGCNA sobre el conjunto de datos GSE54837 utilizando R para identificar módulos relacionados con la EPOC. Antes de la construcción de la red, se identificaron y eliminaron muestras atípicas mediante un análisis de agrupamiento jerárquico utilizando la función hclust con el método de enlace promedio y una métrica de distancia euclidiana. Se seleccionó la potencia óptima de umbralización suave (β = 10) para alcanzar un índice de ajuste a una topología libre de escala R2 ≥ 0,85, equilibrando la topología libre de escala y la conectividad media. Se construyó una matriz de adyacencia y se transformó en una matriz de superposición topológica (TOM). Los módulos génicos se identificaron utilizando el algoritmo dinámico de corte de árboles (deepSplit = 2, minClusterSize = 50). Posteriormente, los módulos con correlaciones de eigengenos > 0,75 se fusionaron utilizando la función mergeCloseModules. Luego, los eigengenos de los módulos se correlacionaron con rasgos clínicos (estado de EPOC, edad, sexo y estado de tabaquismo) mediante coeficientes de correlación de Pearson para identificar módulos asociados a la EPOC destinados a análisis posteriores.
Selección, análisis de enriquecimiento y análisis de la red PPI de genes superpuestos
Se generó un diagrama de Venn utilizando el paquete de R ggvenn para identificar los genes que se superponen entre los DEG, los genes del módulo MEsalmon y los ac4C-RG. Se realizó un análisis de enriquecimiento funcional de los genes superpuestos utilizando las bases de datos Gene Ontology (GO) y Kyoto Encyclopedia of Genes and Genomes (KEGG) con el paquete de R clusterProfiler. Se obtuvo información sobre interacciones proteína-proteína (PPI) a partir de la base de datos STRING (https://string-db.org/) para analizar las interacciones a nivel de proteínas entre los genes superpuestos. El software Cytoscape se utilizó para visualizar la red PPI resultante.
Identificación de genes clave mediante aprendizaje automático
Se aplicaron tres técnicas de aprendizaje automático: regresión con operador de contracción y selección por mínima absoluta (LASSO), potenciación extrema del gradiente (XGBoost) y bosque aleatorio (RF). La regresión LASSO se implementó utilizando el paquete glmnet con validación cruzada de 10 pliegues para determinar el parámetro óptimo de penalización λ. El parámetro type.measure se estableció en «deviance», y el parámetro family se fijó en «binomial». Se seleccionó el valor óptimo de λ mediante el criterio λmin, que minimiza la desviación validada mediante cross-validation, obteniendo 17 genes. XGBoost se realizó utilizando el paquete xgboost con los siguientes hiperparámetros: nrounds = 100, max_depth = 6, eta = 0.3, subsample = 0.8, colsample_bytree = 0.8, y eval_metric = «logloss». La importancia de las características se clasificó según la métrica de ganancia, y se seleccionaron los 30 genes principales. El bosque aleatorio se implementó utilizando el paquete randomForest con ntree = 200. La importancia de las características se ordenó por la disminución media del índice de Gini, y se seleccionaron los 30 genes principales. Los genes seleccionados por los tres métodos de aprendizaje automático se intersectaron para identificar los genes clave para los análisis posteriores.
Construcción y evaluación del modelo de regresión logística para la predicción de riesgo
El conjunto de datos GSE54837 se dividió aleatoriamente en un conjunto de entrenamiento (70 %) y un conjunto de prueba (30 %). Se construyó un modelo de regresión logística sobre el conjunto de entrenamiento utilizando la función glm del paquete MASS, con los niveles de expresión de los genes clave como características de entrada. El rendimiento del modelo se evaluó mediante curvas ROC generadas con el paquete pROC. Los intervalos de confianza del 95 % para el AUC se calcularon mediante 2.000 réplicas bootstrap. La calibración del modelo se evaluó utilizando curvas de calibración generadas con 1.000 remuestras bootstrap (paquete rms). Se realizó un DCA utilizando el paquete dca para evaluar el beneficio clínico neto a lo largo de un rango de probabilidades umbral. Se construyó un nomograma utilizando la función nomogram del paquete rms para facilitar la estimación individualizada del riesgo.
La ecuación de regresión fue:
logit(P) = 0,5823 + 0,6010 × UPP1 - 0,6563 × PTRF + 0,3853 × B4GALT2 - 0,3972 × FAM168B + 0,1848 × PRKCDBP - 0,4787 × TOR3A. (1)
Aquí, P representa la probabilidad predicha de EPOC, y cada coeficiente representa la contribución del valor correspondiente de expresión génica a las probabilidades logarítmicas de EPOC.
Análisis de expresión, red GeneMANIA y red reguladora molecular
Los niveles de expresión génica entre los grupos con EPOC y control en el conjunto de datos GSE54837 se compararon utilizando la prueba de suma de rangos de Wilcoxon. Se generaron gráficos de caja mediante el paquete ggplot2 para visualizar la distribución de los niveles de expresión, con la mediana, el rango intercuartílico (RIC) y los puntos de datos individuales superpuestos. Se utilizó GeneMANIA para construir redes génicas y predecir interacciones funcionales. La búsqueda se realizó con parámetros predeterminados: especie = Homo sapiens, número máximo de genes relacionados = 20. La red resultante se descargó y visualizó, con colores de los enlaces que indicaban los tipos de interacción. Se construyó una red de ARN endógeno competitivo (ceARN) para investigar los mecanismos reguladores postranscripcionales. Las miARN que dianan los seis genes clave se predijeron utilizando dos bases de datos independientes: DIANA-microT (puntuación ≥ 0,8) y miRanda (puntuación ≥ 140, energía ≤ −20 kcal/mol). Se utilizó la intersección de las miARN identificadas por ambas bases de datos para construir pares miARN-mARN. Posteriormente, se predijeron las lncARN que dianan estas miARN mediante la base de datos StarBase. Se construyó y visualizó una red reguladora lncARN-miARN-mARN utilizando Cytoscape. Las relaciones reguladoras transcripcionales se predijeron mediante el análisis de enriquecimiento ChIP-X versión 3 (ChEA3). Para cada gen clave con factores de transcripción predichos, se seleccionaron los 10 principales factores de transcripción con las puntuaciones de enriquecimiento más altas. Se construyó una red reguladora factor de transcripción-diana en Cytoscape.
Análisis de enriquecimiento de conjuntos de genes y evaluación de la infiltración de células inmunitarias
Se realizó un análisis de enriquecimiento de conjuntos de genes (GSEA) utilizando el paquete clusterProfiler para investigar las funciones biológicas de cada gen clave. Para cada gen clave, las muestras se dividieron en grupos de expresión alta y baja según el valor mediano. Se realizó un análisis de expresión diferencial entre los dos grupos utilizando limma, y la lista de genes resultante se ordenó según el cambio en el log₂ con signo. El GSEA se llevó a cabo utilizando la función gseGO para los términos del proceso biológico de GO y la función gseKEGG para las vías de KEGG, con los siguientes parámetros: minGSSize = 10, maxGSSize = 500, pvalueCutoff = 0,05 y nPerm = 1.000. La abundancia relativa de 28 tipos de células inmunitarias se estimó mediante el análisis de enriquecimiento de conjuntos de genes a partir de una única muestra (ssGSEA) implementado en el paquete GSVA. Se obtuvo de estudios previos una matriz de firmas de conjuntos de genes curada que comprende genes marcadores para 28 tipos de células inmunitarias19. Para cada muestra, se aplicó la función gsva con method = "ssgsea", ssgsea.norm = TRUE y kcdf = "Gaussian". Se calcularon los coeficientes de correlación de Spearman entre las puntuaciones de enriquecimiento de ssGSEA y los niveles de expresión de los seis genes clave utilizando la función cor.test. Los valores de p se ajustaron para pruebas múltiples mediante el método de Benjamini-Hochberg. La matriz de correlación se visualizó como un mapa de calor utilizando el paquete pheatmap.
Predicción de fármacos, acoplamiento molecular y análisis de asociación con enfermedades
Se identificaron compuestos terapéuticos potenciales que se dirigen a genes clave utilizando la base de datos DrugBank. Se construyó una red de interacciones «fármaco-dirigido a gen clave» en Cytoscape para visualizar las interacciones predichas entre fármacos y genes. Se realizó el acoplamiento molecular utilizando la plataforma CB-Dock2 para evaluar las afinidades de unión. La estructura tridimensional de la proteína humana UPP1 se obtuvo del Protein Data Bank (PDB ID: 7B8T). Las estructuras moleculares de los fármacos (formato SMILES) se obtuvieron de PubChem. El acoplamiento se realizó utilizando el motor AutoDock Vina, y los resultados se clasificaron según la energía libre de unión (ΔG, en kcal/mol). Los complejos de acoplamiento se visualizaron utilizando PyMOL. Se investigaron las asociaciones entre los genes clave y enfermedades humanas relacionadas con exposiciones ambientales utilizando la base de datos Comparative Toxicogenomics Database (CTD). Cada gen se consultó individualmente, y se extrajeron las diez enfermedades con mayor asociación, las cuales se visualizaron mediante gráficos de radar.
Protocolo de RT-qPCR
Se recogieron muestras de sangre venosa periférica de ocho pacientes con EPOC y ocho controles sanos en el Hospital Shenzhen Luohu de Medicina Tradicional China. El EPOC se diagnosticó según los criterios de la Iniciativa Global para la Enfermedad Pulmonar Obstructiva Crónica (GOLD), definido como una relación VEF1/CVF postbroncodilatador < 0,7. < 0,70. El grupo control estaba compuesto por voluntarios sanos apareados por edad y sexo, sin antecedentes de enfermedades respiratorias y con pruebas de función pulmonar normales (VFE1% predicho ≥ 80% y VFE1/CVF ≥ 0,70). La información basal de los pacientes se muestra en Tabla 2. Se extrajo ARN total de muestras de sangre de EPOC utilizando un kit de extracción de ARN de sangre. Para la síntesis de ADNc, se transcribió inversamente 500 ng de ARN total utilizando un kit de síntesis de ADNc con eliminación de ADN genómico, siguiendo el protocolo proporcionado. El ADNc resultante se diluyó a 150 ng/μL.
Se realizó RT-qPCR utilizando una mezcla maestra de qPCR basada en SYBR Green en un sistema de PCR en tiempo real. Cada reacción de 10 μL contenía 5 μL de mezcla maestra de SYBR Green 2x, 0,5 μL de cada cebador (directo e inverso) (10 μM), 1 μL de cADN diluido (15 ng/μL) y 3 μL de agua libre de nucleasas. Las condiciones de ciclado fueron una desnaturalización inicial a 95 °C durante 5 min, seguida de 40 ciclos de 95 °C durante 10 s y 60 °C durante 30 s, con un análisis final de curva de fusión desde 60 °C hasta 95 °C para verificar la especificidad de la amplificación. Todas las reacciones se realizaron por triplicado técnico. Se utilizó β-actina como gen de referencia interno. La eficiencia de los cebadores para cada gen diana se validó mediante series de dilución de curva estándar y osciló entre el 90 % y el 110 %. Los niveles de expresión génica se normalizaron respecto a β-actina, y la expresión relativa se calculó utilizando el método 2-ΔΔCt. Las comparaciones estadísticas entre los grupos con EPOC y los controles se realizaron utilizando la prueba U de Mann-Whitney.
Análisis estadístico
Las visualizaciones de redes se crearon utilizando Cytoscape, y los análisis estadísticos se realizaron utilizando el software R. A menos que se indique lo contrario, se utilizó la prueba U de Mann-Whitney para datos no normalmente distribuidos, y la prueba t de Student para datos normalmente distribuidos para comparar dos grupos. Un valor de p < 0,05 se consideró estadísticamente significativo.