Critical gene network and signaling pathway analysis of the extracellular signal-regulated kinase (ERK) pathway in ischemic stroke

Microarray datasets related to ischemic stroke were systematically retrieved from the Gene Expression Omnibus (GEO) database1. Three independent microarray expression datasets (GSE22255↗, GSE16561↗, and GSE58294↗) were selected based on stringent inclusion criteria: availability of ischemic stroke samples, matched controls, and high-quality annotation data. Specifically, dataset GSE22255↗ comprised 40 samples, including 20 patients with ischemic stroke and 20 controls, profiled using the Affymetrix Human Genome U133 Plus 2.0 Array (GPL570↗). Dataset GSE16561↗ contained 63 samples (39 ischemic stroke cases and 24 controls), obtained using the Illumina HumanRef-8 v3.0 Expression BeadChip (GPL6883↗). Dataset GSE58294↗ included 92 samples (69 and 23 controls), also profiled on the Affymetrix Human Genome U133 Plus 2.0 Array (GPL570↗).

The ERK signaling pathway genes were comprehensively extracted from the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (KEGG pathway ID: hsa04010). The pathway encompasses a total of 300 annotated genes. These genes were systematically downloaded, and further bioinformatics analyses were conducted to explore their roles and interactions within ischemic stroke pathology. 0002

Three ischemic stroke transcriptomic datasets were retrieved from the GEO. Raw RNA expression data and platform annotation files (GPL) were obtained using the GEOquery package in R (v2.68.0). Gene symbols were annotated for each dataset separately. To ensure a robust and integrative analysis, datasets GSE22255↗ and GSE16561↗ were merged into a combined discovery cohort, while GSE58294↗ served as an independent validation set. Batch effects introduced by different sequencing platforms were corrected using the ComBat function from the sva package in R.

To evaluate the distribution of gene expression before and after batch effect correction, box plots of the merged dataset were generated. Additionally, uniform manifold approximation and projection (UMAP) analysis was performed using the umap package in R to visualize the overall structure of the dataset before and after batch effect removal. This analysis facilitated the assessment of clustering patterns among ischemic stroke and control samples, ensuring the integrity and comparability of the integrated dataset.

Following batch effect correction, a combined dataset was obtained consisting of 103 samples, including 44 healthy controls and 59 ischemic stroke samples. To facilitate downstream comparisons, all healthy controls were placed at the beginning of the dataset, followed by the diseased samples. The analysis of differential expression genes was performed using the limma package in R. Genes were considered differentially expressed if they met the following criteria: an adjusted p-value (Benjamini–Hochberg correction) of <0.05 and an absolute log2 fold change (|log2FC|) of >0.3. This threshold corresponds to an approximately 1.23-fold change in expression level.

A volcano plot was generated using the ggplot2 package to depict the statistical significance (−log10 of the adjusted p-value) against the magnitude of expression change (log2 fold change). Vertical dashed lines were drawn at −0.3 and 0.3 on the x-axis, and a horizontal dashed line at −log10 (0.05) on the y-axis, demarcating the thresholds for significance. Additionally, significantly up- and downregulated genes were subsequently visualized using a heatmap generated with the pheatmap package.

A WGCNA was performed using the set of differentially expressed genes (DEGs) derived from the integrated, batch-corrected dataset (n = 103 samples). First, hierarchical clustering was conducted using average linkage to detect outlier samples. Four samples—GSM416554↗, GSM416550↗, GSM416535↗, and GSM416539↗—were excluded from subsequent analyses based on their extreme dissimilarity to the main cluster. This filtering step ensured a more coherent and reliable input dataset for network construction.

Next, the pickSoftThreshold function in the WGCNA package was used to identify an optimal soft thresholding power for the network. The power value of 6 was selected, as it satisfied the approximate scale-free topology criterion and yielded a high mean connectivity that balanced sensitivity and specificity in module detection.

Using a power of 6, the network was built by calculating pairwise correlations between the expression profiles of the DEGs and transforming them into a weighted adjacency matrix. The minimum module size (minModuleSize) was set to 40 to exclude small, potentially spurious modules, and modules with similar eigengenes were merged using a mergeCutHeight of 0.25. Each module was assigned a unique color identifier for clarity. Ultimately, this clustering procedure yielded two main co-expression modules.

Module–trait relationships were evaluated by correlating each module eigengene with the sample classification (healthy vs. ischemic stroke). Modules that did not correlate significantly with the trait were grouped into the “gray” category and excluded from further consideration. Genes within the modules of interest were assessed for both module membership (MM) and their involvement in the ERK pathway. The ggvenn package was employed to intersect module-member genes with ERK-related genes, identifying the most relevant candidates—termed “key genes”—for subsequent analyses.

Gene Ontology (GO) enrichment analysis was conducted to characterize the biological roles of the gene signature ERK (GSERK) genes. The clusterProfiler package and org. Hs.eg.db in R were used to calculate enrichment scores across three GO domains: biological process (BP), cellular component (CC), and molecular function (MF). For the BP and MF categories, genes were considered significantly enriched if they met thresholds of p-value < 0.05 and q-value cut off < 0.05. Because GSERK genes showed relatively weaker enrichment in the CC domain, we relaxed the criteria to p-value < 0.1 and q-value cut off < 0.1 for CC enrichment analyses. Enrichment results were visualized using bubble plots and circle plots to highlight the most prominently enriched GO terms.

The identified GSERK genes were uploaded to the STRING databaseto predict protein–protein interactions (PPI). The resulting network was exported and further visualized in Cytoscape, where edges represent putative functional or physical associations among the GSERK proteins. 0003

To prioritize key genes within the GSERK network, the cytoHubba plug-in in Cytoscape was employed. The computational algorithms were applied sequentially to rank the genes by their topological properties. Any gene not ranked as a hub by any of these methods was excluded from final consideration.

The hub GSERK genes were selected for in-depth expression profiling in ischemic stroke relative to healthy control samples. Expression data for these genes were extracted from both the complete dataset (“total sample set”) and an independent validation cohort. Boxplots were generated with ggpubr and forcats in R language, stratifying samples into ischemic stroke versus healthy control groups. Statistical significance between the two groups was determined using appropriate post hoc tests (e.g., t-tests or non-parametric equivalents), with p-values adjusted for multiple comparisons where necessary.

To evaluate the potential diagnostic utility of each hub GSERK gene, ROC curves were plotted and area under the curve (AUC) values were calculated using the pROC package. ROC analysis was performed separately for the total sample set and the validation cohort, with higher AUC values indicating stronger discriminatory power between IS and healthy subjects.

To robustly discern key biomarkers for ischemic stroke, four complementary machine learning algorithms were applied to the gene expression dataset. First, the Boruta algorithm was used to iteratively assess feature relevance by comparing actual predictor importance scores to those of randomized features. Second, a support vector machine (SVM) approach was used with a recursive feature elimination strategy, wherein features were iteratively removed based on their contribution to classification accuracy. Third, a least absolute shrinkage and selection operator (LASSO) model was implemented to impose an L1 penalty on regression coefficients, effectively shrinking the influence of less informative genes and retaining only the most predictive variables. Finally, a random forest classifier evaluated feature importance via the mean decrease in the Gini index, prioritizing features with the highest impact on classification purity at each node.

To maximize predictive power and minimize algorithm-specific biases, we compared the feature sets identified by each of the four methods. The final gene list was derived by taking the intersection of all selected features, thereby pinpointing genes that consistently emerged as top predictors across multiple machine learning paradigms.

To delineate immune heterogeneity in ischemic stroke, we stratified the harmonized dataset (n = 103) into two consensus clusters (k = 2) using K-means clustering with 1,000 iterations and Euclidean distance. Cluster stability was validated via consensus matrix analysis and cumulative distribution function (CDF) curves. Immune cell infiltration was quantified via single-sample gene set enrichment analysis (ssGSEA) using 28 immune cell type-specific gene sets. Stromal and immune scores were computed using the ESTIMATE algorithm, and inter-cluster differences were assessed via Wilcoxon rank-sum tests.

To establish a clinically translatable risk assessment framework for ischemic stroke, we constructed a nomogram integrating the GSERK hub genes using multivariable logistic regression. Gene expression values were standardized (z-scores) and assigned weighted points proportional to their regression coefficients. The model was validated via bootstrapping (1,000 resamples) to estimate calibration and discrimination metrics. Odds ratios (ORs) and 95% confidence intervals (CIs) for each gene were derived from univariate and multivariate analyses. Statistical significance was assessed using Wald tests, with p < 0.05 considered significant.

To evaluate the expression levels of ERK pathway–related proteins DUSP1 and GADD45A under hypoxia/reoxygenation (H/R) conditions, SH-SY5Y neuroblastoma cells were cultured and subjected to the following experimental treatments: (1) control group: Cells were maintained under normoxic conditions (21% O₂, 5% CO₂) for the full duration; (2) H/R group: cells were exposed to hypoxia (1% O₂, 5% CO₂, 94% N₂) for 6 h, followed by reoxygenation under normoxic conditions (21% O₂) for 24 h.

Following treatment, total protein was extracted using RIPA buffer supplemented with protease and phosphatase inhibitors. Protein concentrations were quantified using the BCA protein assay kit (Thermo Fisher Scientific, United States). Equal amounts of protein (20–30 μg) were separated on SDS-PAGE gels and transferred onto PVDF membranes. Membranes were blocked in 5% non-fat milk and incubated overnight at 4°C with the following primary antibodies: anti-DUSP1 (1:1000, Thermo Fisher Scientific, TA890036), anti-GADD45A (1:1000, Abcam, ab180768), and anti-GAPDH (1:5000, Thermo Fisher Scientific, MA5-15738) as a loading control. After washing, membranes were incubated with HRP-conjugated secondary antibodies (1:5000) for 1 h at room temperature. Protein bands were visualized using ECL detection reagents and quantified by ImageJ software. Relative expression levels of DUSP1 and GADD45A were normalized to GAPDH.

To evaluate the effectiveness of batch effect correction across two ischemic stroke transcriptomic datasets (GSE16561↗ and GSE22255↗), we assessed their expression distribution patterns before and after data harmonization. Principal component analysis (PCA) revealed strong dataset-specific clustering in the uncorrected data, along with divergent median expression values (Figure 1A). After applying a robust batch correction algorithm, these distributions converged substantially, and median values aligned more closely, indicating effective mitigation of inter-dataset variability (Figure 1B).

We further validated these findings using UMAP. In the uncorrected data, GSE16561↗ and GSE22255↗ occupied distinct, non-overlapping clusters, reflecting pronounced technical differences between the cohorts (Figure 1C). Post-correction, samples from both datasets became interspersed, suggesting that batch-driven artifacts were minimized while preserving underlying biological variation (Figure 1D). These observations underscore the importance of implementing rigorous batch correction strategies for integrative analyses of heterogeneous genomic datasets in ischemic stroke research.

From the total of 103 samples subjected to differential expression analysis (44 healthy and 59 diseased), 140 genes met the significance criteria (adjusted p-value <0.05 and |log2FC| > 0.3). Of these, 125 genes were upregulated in ischemic stroke samples compared with healthy controls, whereas 15 genes were downregulated. The volcano plot revealed that distribution of DEGs in Figure 2A. Additionally, heatmap analysis of the DEGs further illustrated distinct expression profiles between healthy controls and ischemic stroke samples (Figure 2B).

Hierarchical clustering revealed four outlier samples (GSM416554↗, GSM416550↗, GSM416535↗, GSM416539↗), which were removed to ensure the resulting network captured genuine co-expression relationships rather than artifacts from disparate data points (Figure 3A). By systematically evaluating the scale-free topology fit at different power values, a soft threshold of 6 was chosen (Figure 3B). Under this optimal threshold, co-expression modules were identified using dynamic tree cutting with minModuleSize set to 40. Subsequent merging of closely related modules (mergeCutHeight = 0.25) resulted in two major co-expression modules, each assigned a distinct color tag (Figure 3C). Correlation analysis between module eigengenes and sample phenotypes indicated that at least one of the resulting modules exhibited strong relevance to the ischemic stroke phenotype. A heatmap of module–trait correlations confirmed this module as a primary focus of further investigation (Figure 3D).

A scatter plot of module membership versus gene significance (or correlation to the trait) indicated a predominantly negative correlation for genes in the key module (Figure 3E), suggesting potential modulatory or inhibitory roles in the disease process. Using the ggvenn package, six ERK-related genes (GADD45A, DUSP1, EREG, IL1beta, JUN, and GADD45B) were found at the intersection of module membership and known ERK pathway components (Figure 3F). These genes, collectively referred to as “gene signature ERK pathway (GSERK),” emerged as critical hubs for further functional validation, laying the groundwork for subsequent exploration of their roles in ischemic stroke pathophysiology.

Under the specified thresholds (p < 0.05, q < 0.05 for BP and MF), the GSERK gene set displayed significant enrichment in 163 biological processes and 13 molecular functions. For the cellular component category, nine enriched terms were identified using a more permissive cutoff (p < 0.1, q < 0.1).

In biological process enrichment, a total of 32 BP terms passed a more stringent cutoff of p < 0.01 and q < 0.01 (Figures 4A,B). Notably, GSERK demonstrated the strongest enrichment in regulation of p38MAPK cascade, p38MAPK cascade itself, regulation of protein kinase activity, and regulation of kinase activity. Each of these processes involved five of the six GSERK members (GADD45A, DUSP1, IL1B, GADD45B, and EREG). These findings underscore a mechanistic link between GSERK and MAPK-driven signaling pathways, which are known to be pivotal in ischemic stroke pathophysiology.

While GSERK genes did not exhibit strong enrichment in the CC category under the standard significance thresholds, we identified nine enriched cellular components at p < 0.1 and q < 0.1 (Figures 4C,D). This less stringent criterion allowed the capture of potential CC associations that may still contribute to the spatial dynamics of these critical genes within the cell.

Among the 13 significantly enriched MF terms (Figures 4E,F), growth factor receptor binding pathway showed the highest enrichment ratio under MF enrichment. EREG and IL1B were specifically implicated in this functional category, suggesting they may modulate key receptor-mediated signaling events pertinent to ERK pathway activation in ischemic stroke. Collectively, these GO analyses point to a central role for GSERK genes in modulating kinase-related signaling cascades and growth factor receptor interactions, further highlighting their potential as therapeutic targets or biomarkers in ischemic stroke research.

The STRING-based PPI network (Figure 5A) highlighted direct and indirect links among the six GSERK proteins, indicating extensive interconnectivity within stress- and inflammation-related pathways. Cytoscape-based refinement (Figure 5B) enabled more detailed visualization, revealing strong central positioning for five of the GSERK genes (GADD45A, DUSP1, IL1B, JUN, GADD45B), whereas EREG showed fewer high-confidence interactions.

Applying seven distinct cytoHubba algorithms, EREG did not meet the criteria for a hub gene in any of the methods employed. Consequently, the remaining five genes were designated as the hub GSERK members: GADD45A, DUSP1, IL1B, JUN, and GADD45B. Color-coded network diagrams for each algorithm (Figures 5C–I) consistently highlighted these five genes as central nodes, underscoring their potential significance in modulating ERK-related signaling cascades implicated in ischemic stroke pathogenesis.

In the total sample set, all five hub GSERK genes were upregulated in IS samples compared with healthy controls (Figure 6A). Of particular note, GADD45A, DUSP1, and GADD45B showed marked differences (p < 0.001) between the two groups. In the validation cohort (Figure 6C), four genes (GADD45A, DUSP1, GADD45B, and IL1B) were significantly upregulated in IS relative to healthy controls (p < 0.001), whereas JUN did not display a statistically significant difference. These findings highlight potential roles for GADD45A, DUSP1, GADD45B, and IL1B as robust markers of ischemic stroke.

The ROC curve analysis in the total sample set revealed moderate discriminatory capacity for DUSP1 (AUC = 0.75), GADD45A (AUC = 0.72), and GADD45B (AUC = 0.70), suggesting these genes may serve as potential diagnostic biomarkers for ischemic stroke (Figure 6B). Consistent with these results, the validation cohort showed high diagnostic performance for DUSP1 (AUC = 0.91), followed by GADD45A (AUC = 0.79) and GADD45B (AUC = 0.72) (Figure 6D). Collectively, these results support the potential of DUSP1, GADD45A, and GADD45B as useful biomarkers for ischemic stroke, warranting further clinical validation.

All four machine learning algorithms identified a partially overlapping set of candidate biomarkers, reflecting each method’s distinct way of handling noise and collinearity in the data (Figures 7A–E). Upon intersecting the results (Figure 7F), three genes—GADD45A, DUSP1, and GADD45B—were consistently ranked as important across Boruta, SVM, LASSO, and random forest analyses. These robust markers thus hold particular promise for further validation and potential translational applications in ischemic stroke research.

To resolve immune-stromal heterogeneity in ischemic stroke, we stratified the harmonized cohort (n = 103) into two molecular subgroups using K-means clustering (k = 2, Figure 8A). Consensus matrix analysis confirmed robust cluster stability, as indicated by an area under the CDF curve (Figure 8B). Subgroup 2 exhibited a pronounced hyperinflammatory phenotype, characterized by elevated infiltration of effector memory CD8⁺ T cells, activated CD8⁺ T cells, macrophages, and neutrophils. In contrast, Subgroup 1 showed enrichment for natural killer cells, T helper 17 (Th17) cells, and central memory CD8⁺ T cells (Figure 8C). These findings suggest divergent immune microenvironments that could influence therapeutic responsiveness.

Application of the ESTIMATE algorithm further distinguished the two subgroups. Subgroup 1 demonstrated significantly higher StromalScores (p = 0.009) and ESTIMATEScores (p = 0.009), indicative of active immune-stromal crosstalk (Figure 8D). In contrast, Subgroup 1 also showed muted stromal remodeling (p = 0.4), potentially reflecting adaptive mechanisms to aiming at mitigating ischemic damage. Notably, these subgroup-specific profiles align with ERK pathway dynamics, wherein hyperinflammation may exacerbate oxidative stress, whereas more quiescent stromal remodeling could temper reparative signaling.

To translate GSERK gene signatures into a clinically applicable model, we constructed a nomogram incorporating the expression levels of GADD45A, DUSP1, and GADD45B, along with age and gender as clinical covariates (Figure 9A). Each variable contributed a weighted score based on multivariable logistic regression coefficients, and the total score was mapped to an estimated stroke risk probability (ranging from 0.1 to 0.9). The nomogram offers a visual, interpretable tool for individualized risk assessment. Multivariate logistic regression analysis (Figure 9B) revealed that among the molecular variables, DUSP1 (OR = 2.39, 95% CI: 0.87–6.57; p = 0.092) and GADD45B (OR = 2.53, 95% CI: 0.63–10.22; p = 0.192) showed a trend toward stroke association, though neither reached statistical significance. GADD45A exhibited a weaker effect (OR = 1.92, 95% CI: 0.70–5.71; p = 0.298). Notably, age emerged as a statistically significant predictor of stroke risk (OR = 1.05, 95% CI: 1.01–1.08; p = 0.015), highlighting the importance of integrating clinical parameters into molecular models. Gender was not significantly associated with outcome (male vs. female: OR = 1.05, p = 0.922). These results suggest that while ERK-associated genes, particularly DUSP1 and GADD45B, may contribute to ischemic stroke susceptibility, age remains the dominant independent predictor in this integrative model. Further validation in larger cohorts may assist clarify the prognostic relevance of GSERK components.

Western blot analysis revealed that both DUSP1 and GADD45A protein levels were significantly upregulated in SH-SY5Y cells subjected to H/R compared to normoxic controls (Figure 10A). Densitometric quantification demonstrated an increase in DUSP1 expression (p < 0.01, Figure 10B) and GADD45A expression (p < 0.05, Figure 10B), normalized to GAPDH. These findings support the transcriptomic prediction that ERK-associated genes are responsive to ischemic-like stress and may play a regulatory role in the cellular response to oxidative injury.