Identification and validation of glucocorticoid receptor and programmed cell death-related genes in spinal cord injury using machine learning

Abstract

Spinal cord injury (SCI) is a severe neurological disorder, with glucocorticoids like methylprednisolone commonly used for treatment. However, their efficacy and risks remain controversial. Programmed cell death (PCD) mechanisms have been increasingly implicated in SCI pathology. This study aimed to identify differentially expressed genes (DEGs) related to glucocorticoid receptors and PCD and to construct a diagnostic model to guide glucocorticoid use in SCI treatment. SCI datasets (GSE5296, GSE47681, GSE151371, and GSE45550) were analyzed using protein-protein interaction networks, consensus clustering, GSVA for PCD pathway enrichment, and WGCNA. A total of 113 diagnostic models were developed through 12 machine learning algorithms, with the optimal model, "Lasso + Stepglm[both]," featuring six genes: Abca1, Cdh1, Glipr1, Glt8d2, Il10ra, and Pde5a. Validation through qRT-PCR confirmed the differential expression of four genes (Abca1, Glipr1, Il10ra, and Cdh1), which demonstrated strong predictive performance. Pathway enrichment of GRRDEGs was analyzed using GO, KEGG, and Bayesian network methods, and immune cell infiltration was assessed via CIBERSORT. In this study, we identified GR- and PCD-related DEGs in SCI and constructed a diagnostic model that may improve understanding of SCI molecular mechanisms and inform future investigations of glucocorticoid use.

Fig. 1

Technology roadmap. DEGs, Differentially Expressed Genes; GRRGs, Glucocorticoid receptor-Related Genes; GRRDEGs, Glucocorticoid receptor-Related Differentially Expressed Genes; GO, Gene Ontology; GSVA, Gene Set Variation Analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; ROC, Receiver Operating Characteristic; SCI, spinal cord injury; PPI, Protein-protein Interaction; WGCNA, Weighted Gene Co-expression Network Analysis;

Fig. 2

Differential gene expression analysis. (A) Volcano plot illustrating the differential expression analysis of SCI and Control samples in the combined GEO datasets. (B) Volcano plot showing the analysis of GRRDEGs in the integrated GEO dataset. (C) Heat map displaying GRRDEGs in the integrated GEO datasets.

Fig. 3

PPI interaction network and GO/KEGG enrichment analysis. (A) The PPI network of GRRDEGs calculated using the STRING database. (B) Venn diagram of the top 20 GRRDEGs from five algorithms in the CytoHubba plugin. (C–D). Bar plot (C) and bubble plot (D) showing the results of GO and KEGG enrichment analysis of GRRDEGs. Pathway diagram adapted from KEGG (Kyoto Encyclopedia of Genes and Genomes)^–.

Fig. 4

Identification of GRRDEGs in spinal cord injury. (A) Box plot showing the expression profiles of 17 GRRDEGs. (B) Heatmap depicting the expression patterns of 17 GRRDEGs between SCI and Control samples. (C) Chromosomal locations of the 17 GRRDEGs. (D) Gene relationship circle diagram for the 17 GRRDEGs, with violet lines indicating positive correlations and green lines indicating negative correlations. (E) Relationships among the 17 GRRDEGs. Significance levels were indicated by *p < 0.05, **p < 0.01, ***p < 0.001.

Fig. 5

Combined datasets immune infiltration analysis by CIBERSORT algorithm. (A–B) The proportion of immune cells in theCombined Datasets bar graph (A) and group comparison graph (B). (C) Correlation heatmap showing immune cell infiltration abundance in the combined datasets. (D) Bubble plot illustrating the correlation between GRRDEGs and immune cell infiltration abundance in the combined datasets. * represents p value < 0.05, statistically significant; ** represents p value < 0.01, highly statistically significant; *** represents p value < 0.001 and highly statistically significant. In the proportion bar chart and group comparison chart, violet is SCI samples, and green is Control samples. green is negative correlation, violet is positive correlation, and the depth of color represents the strength of correlation.

Fig. 6

Consensus clustering based on GRRDEGs expression matrices. (A) Consensus clustering matrix for k = 2. (B) Relative change in the area under the CDF delta curves. (C) Cumulative distribution function (CDF) for consensus clustering at k = 2–6. (D) Consensus clustering scores. (E) PCA plot of the two clusters, with each scatterplot point representing a sample based on the two principal components (PC1 and PC2) derived from GRRDEG expression profiles.

Fig. 7

Gene abundances and immune characteristics of GRRDEGs. (A) and (B) display the expression patterns of 17 GRRDEGs across two distinct glucocorticoid receptor-related phenotypes using boxplots and heatmaps, respectively. (C) Boxplot showing the proportions of 22 immune-infiltrated cell types analyzed via the CIBERSORT algorithm. (D) Bubble plot illustrating the correlation between GRRDEGs and immune cell infiltration abundance. Significance levels were indicated by *p < 0.05, **p < 0.01, ***p < 0.001.

Fig. 8

Bayesian network enrichment analysis. Bar plot (A) and bubble plot (B) showing the results of Bayesian network enrichment analysis of GRRDEGs. (C) Bayesian network inference of the pathway regulatory network of GRRDEGs. Each node represents a pathway; the redder and larger the node, the higher the gene expression value. The thicker the line, the stronger the regulatory effect.

Fig. 9

GSVA analysis of PCD pathways and Identification of the critical gene modules. (A) GSVA analysis of 13 PCD pathways displayed as box plots comparing Control and SCI groups. (B) Module-trait correlations, with columns representing traits and rows representing module eigengenes. (C) The correlation between cluter and module eigengenes. Each column matches a trait, and each row corresponds a module eigengene. Each unit cell includes the corresponding p-value and correlation. SCI, Spinal Cord Injury.

Fig. 10

Diagnostic performance of our model. (A) Volcano plot of the included genes. (B) ROC curves for the included genes. (C) qRT-qPCR validation of GRPRDEG mRNA expression. (D) Nomogram predicting the risk of SCI based on GRPRDEGs (Abca1, Glipr1, Il10ra). (E–G) ROC curves evaluating the prediction model in the combined dataset, GSE151371 dataset, and GSE45550 dataset.