The gene expression datasets GSE124272 and GSE150408 were obtained from the Gene Expression Omnibus (GEO) hosted by the National Center for Biotechnology Information (https://www.ncbi.nlm.nih.gov/geo/) as detailed in Table 1. The original GEO records defined the patient selection criteria for datasets GSE150408 and GSE124272 as follows. The sciatica groups comprised individuals diagnosed with either sciatica or lumbar disc prolapse, both confirmed by magnetic resonance imaging (MRI). The healthy control groups consisted of volunteers without low back pain or sciatica. To maximize data utility for model development while ensuring an unbiased performance estimate, we employed a pre-defined, dataset-level analytical strategy: data from both GSE124272 and GSE150408 were integrated for model development and feature selection, while the entire GSE150408 cohort was reserved as a fixed, hold-out internal validation set for final model evaluation. The gene expression profiling in this study was based on data generated using the Agilent-072363 SurePrint G3 Human GE v3 8 × 60K Microarray platform (GPL21185). The GSE124272 dataset was extracted and divided into 8 healthy normal individuals and 8 individuals with sciatica. For the GSE150408 dataset, we specifically used the baseline (pre-treatment) samples, comprising 17 cases from the normal group and 17 from the sciatica group, prior to any therapeutic intervention. The GSE150408 dataset provided 17 cases from the normal group and 17 from the sciatica group. The pre-processed “Series Matrix” files (e.g., GSE124272_series_matrix.txt) were downloaded from GEO. These files contain background-corrected and summarized expression values. After integrating the datasets from GSE124272 and GSE150408, probes were annotated using the platform-specific annotation file (GPL21185). Quantile normalization was then performed using the normalizeBetweenArrays function from the limma package (version 3.58.1) to ensure comparability across all samples (Ritchie et al., 2015). Subsequently, to adjust for technical batch effects between the two studies, the ComBat function from the sva package (version 3.46.0) was applied to the normalized data using the parametric empirical Bayesian framework (par.prior = TRUE) (Leek et al., 2012). Principal-component analysis was performed on the merged expression matrix before and after ComBat; samples were visualised by dataset to inspect clustering patterns. The proportion of variance attributable to batch was quantified as the ratio of the between-dataset sum of squares to the total sum of squares on the first two PCs (PC1 + PC2). The resulting batch-corrected expression matrix was used for all downstream analyses. Genes associated with neuroinflammatory mechanisms were identified by querying the term “neuroinflammation” across all annotation sections (“All Sections”) and restricting to the “Protein Coding” category in the GeneCards database (https://www.genecar-ds.org/, accessed on May 17, 2025). An initial broad set of genes was obtained from this query, to which a stringent relevance score filter (>1.5) was applied, yielding the final gene set used for subsequent analysis.

We screened for differentially expressed genes (DEGs) that were significantly expressed in both the control and sciatica groups. Genes with logFC >0.58 and p < 0.05 are considered upregulated genes. Conversely, genes with logFC <−0.58 and p < 0.05 are classified as downregulated. The chosen criteria for gene expression analysis are |log2(FC)| >0.58 and p < 0.05, as shown in a volcano plot. Subsequently, we employed the pheatmap package (version 1.0.13) to illustrate and evaluate gene expression differences between the normal and the sciatica groups.

The analysis of Gene Ontology (GO) functional annotations for DEGs related to neuroinflammation was performed using the “clusterProfiler” package (version 4.0) in the R programming environment (Yu et al., 2012). The gene set ‘c2. cp.v7.5.1. symbol-s. gmt’ was sourced from the Molecular Signatures Database (MSigDB, version 2025.1. Hs, accessed on May 17, 2025) and subsequently employed to conduct Gene Set Enrichment Analysis (GSEA) on the identified DEGs (Liberzon et al., 2015).

Machine learning analysis followed a pre-defined, two-stage workflow (development/tuning and final evaluation) to ensure robustness. Prior to model training, all feature variables were centered and scaled to a zero mean and unit variance. In the first stage (model development and feature selection), expression data from the integrated GSE124272 and GSE150408 cohorts were pooled. Model development and tuning were conducted using k-fold cross-validation on this combined dataset. For least absolute shrinkage and selection operator (LASSO) regression, performed using the “glmnet” package (version 4.1) (Friedman et al., 2010), the optimal regularization parameter (λ) was determined through 10-fold cross-validation with a fixed random seed (seed = 2022). For feature selection, we employed the Support Vector Machine-Recursive Feature Elimination (SVM-RFE) methodology (Duan et al., 2005) with a radial basis function kernel; its hyperparameters (cost = 0.5, gamma = 0.015625) were optimized via 5-fold cross-validation with a fixed random seed (seed = 2024). The overlapping key genes identified by both methods were extracted using the “Venn” package. In the second stage (final model evaluation), the diagnostic model was constructed using only the consensus hub genes identified above. This final model was then applied, for the first and only time, to the pre-specified, hold-out GSE150408 validation cohort to generate unbiased performance metrics. Calibration curve analysis was performed using the rms package (version 6.4-0) and ResourceSelection package (version 0.3-5). Decision curve analysis was performed utilizing the rmda package (version 1.6) (Vickers and Elkin, 2006). The area under the receiver operating characteristic curve (AUC) and its 95% confidence interval (CI) for the model (Figure 7) were calculated using the pROC R package (version 1.18.0), which implements DeLong’s method. All statistical analyses were performed using R version 4.2.1.

CIBERSORTx employs linear support vector regression to decipher transcriptomic expression matrices, enabling the estimation of the composition and relative abundance of immune cell types within a sample mixture. The gene expression matrix data were analyzed using the standard CIBERSORT algorithm with the LM22 signature matrix (547 genes defining 22 immune cell types) and the default setting of 1 permutation for significance analysis (Newman et al., 2019). This method estimates immune cell fractions, retaining only samples with a CIBERSORT inference p-value <0.05 to ensure reliable deconvolution results, thereby generating a final immune cell infiltration matrix. The R software package “ggplot2” (version 3.4.2) was used to produce hist-ograms illustrating the distribution of immune cell types within each sample. To evaluate the association between immune cell infiltration and DEGs, the Spearman correlation coefficient used to generate a correlation heatmap. Additionally, to explore disparities in immune cell infiltration across groups, boxplots were constructed, aiding the identification of immune infiltration patterns associated with sciatica. To ensure robustness, we performed cross-validation using the MCP-counter algorithm.

We retrieved miRNAs linked to key genes from the TargetScanHuman database (https://www.targetscan.org/vert_80/) and subsequently constructed a regulatory network depicting interactions between mRNAs and miRNAs. We used the ChIPBase v3.0 database (https://rnasysu.com/chipbase3/index.php) to identify transcription factors that bind to hub genes. We constructed an interaction network linking key genes to transcription factors (TFs) after retrieving relationships between TF targets and DEGs. To identify potentially repurposable drugs, we queried DGIdb (version 5.0) using our four hub genes (PLG, LRRK2, NLRP3, KLRK1). We applied the “Approved” drug status filter available in the database interface to prioritize clinically translatable agents. Cytoscape software (version 3.9.1) was then used to visualize the mRNA-miRNA regulatory network along with the mRNA-TF and mRNA-drug interaction networks.

Peripheral blood samples were obtained from 15 patients with MRI-confirmed lumbar disc herniation presenting radiating leg pain persisting >3 months and scheduled for lumbar discectomy (8 male/7 female, age 60.2 ± 6.8 years), along with 15 age- and sex-matched healthy controls (age 63.9 ± 7.2 years) without history of chronic pain or neurological disorders. Total RNA was extracted from ≤200 μL whole blood using the EasyPure® miRNA Kit (TransGen Biotech). RNA concentration and purity (A260/280) were measured with a NanoVue Plus spectrophotometer (Cytiva). Reverse transcription was performed with 1 μg total RNA using the Evo M-MLV RT Mix Kit with gDNA Clean for qPCR Ver.2 (Accurate Biology) in a 20 μL reaction (37 °C, 15 min; 85 °C, 5 s). qPCR was run on a QuantStudio 1 system (Thermo Fisher) with PerfectStart® Green qPCR SuperMix (TransGen Biotech). Each 20 μL reaction contained 10 μL SuperMix, 2 μL cDNA (40× diluted), 0.4 μM each primer, and nuclease-free water. Cycling conditions: 94 °C for 2 min; 45 cycles of 94 °C for 5 s and 60 °C for 30 s; melt-curve from 65 °C to 95 °C (0.5 °C/5 s). Gene expression was normalized to ACTIN and calculated using the ΔΔCt method. All primers (designed with Primer-BLAST, amplicons 80–150 bp, Tm ≈ 60 °C) are listed in Supplementary Table S1.

Data calculations and analyses were performed using R (version 4.3.3). The normality of continuous variables was assessed using the Shapiro-Wilk test, and the homogeneity of variances was evaluated using Levene’s test. Based on these results, comparisons between two groups for continuous variables were conducted with the Student’s t-test or the Mann-Whitney U test, with p-values <0.05 considered statistically significant. To control for multiple hypothesis testing in high-throughput analyses, the Benjamini–Hochberg false discovery rate (FDR) correction was applied to the p-values from the genome-wide differential expression analysis. In contrast, for the Gene Ontology (GO) functional enrichment analysis and Gene Set Enrichment Analysis (GSEA), the algorithm-generated FDR-corrected q-values are reported, with the complete FDR-adjusted results provided in the supplementary materials for full transparency (Supplementary Table S2). Data from qRT-PCR are presented as mean ± standard deviation (SD). Differences in gene expression levels between the sciatica patient group and the healthy control group were assessed for statistical significance using the two-tailed Student’s t-test. P < 0.05 was considered statistically significant. The overall study design is summarized in Figure 1.

After processing the dataset, we obtained 25 samples from the healthy normal group and 25 samples from the sciatica group. The distributions of the data were plotted before and after normalization (Figures 2A,B). Before batch adjustment, the first two principal components (PC1 = 28.1%, PC2 = 9.8%) displayed pronounced dataset-specif- ic segregation that explained 16.2% of the total variance (Supplementary Figure S1; Supplementary Table S3), indicating substantial technical bias. After ComBat correction, this batch-associated variance fell to 2.1% and the dataset clusters vanished (Figure 2C; Supplementary Table S4). We identified 449 neuroinflammation-related genes from GeneCards (https://www.genecards.org/). The combined MSigDB neuroinflammation keyword search resulted in a reference set of 185 unique genes. Notably, 94% of these genes (174 out of 185) were also present in our GeneCards-derived set (Supplementary Table S5). This striking overlap strongly supports the comprehensiveness and specific biological relevance of our initial gene selection for studying neuroinflammation. Principal Component Analysis of the samples from both groups revealed distinct variations in the distribution patterns of 351 associated genes. The differentially expressed genes were shown using volcano plots, with 206 genes upregulated and 145 genes downregulated (Figure 2D). We identified 11 DEGs related to neuroinflammation that differed significantly between sciatica patients and healthy controls (Figure 3A). The chromosomal locations of the eleven DEGs are shown in Figure 3B. In contrast, Figure 3C presents a heatmap of neuroinflammation gene expression differences between the sciatica and normal groups, and Figure 3D shows the correlation among these genes using Spearman correlation analysis.

GO analysis of differentially expressed neuroinflammation genes in sciatica is shown in Table 2. The GO analysis indicated enrichment in various biological processes, including: inflammatory response, immune effector processes, leukocyte migration, cytokine production, and lymphocyte-mediated immunity (Figures 4A–C). Additionally, it highlighted the following cellular components: the inflammasome complex, the specific granule lumen, the external side of the plasma membrane, the autolysosome, and the secretory granule lumen (Figures 4A,B,D). It also identified various molecular functions: immune receptor activity, pattern recognition receptor activity, MAP kinase kinase kinase activity, ionotropic glutamate receptor activity, and peptidoglycan binding (Figures 4A,B,E). Table 3 presents the GSEA results related to sciatica. GSEA results showed positive enrichment of genes in the REACTOME neutrophil degranulation pathway (Normalized Enrichment Score (NES) = 2.638, FDR <0.001) (Figure 5A) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) complement and coagulation cascades pathway (NES = 1.915, FDR <0.001) (Figure 5B). REACTOME interferon gamma signaling (NES = 2.095, FDR <0.001) (Figure 5C), KEGG toll-like receptor signaling pathway (NES = 1.770, FDR <0.001) (Figure 5D), and KEGG cytokine-cytokine receptor interaction (NES = 1.561, FDR <0.001) (Figure 5E). In contrast, negative enrichment was observed in several pathways, including REACTOME DNA repair (NES = −1.629, FDR <0.001) (Figure 5F), KEGG cell cycle (NES = −1.731, FDR <0.001) (Figure 5G), and KEGG ribosome (NES = −2.510, FDR <0.001) (Figure 5H), REACTOME translation (NES = −2.633, FDR <0.001) (Figure 5I). The findings indicate that the expression levels of various neuroinflammation genes are significantly associated with both the initiation and advancement of sciatica.

We analyzed the DEGs in the dataset using two algorithms. The LASSO regression identified six key genes (Figures 6A,B) (Supplementary Tables S6, S7), with the optimal regularization parameter (λ = 0.082924) selected through 10-fold cross-validation minimizing binomial deviance (Supplementary Table S8). The SVM machine learning algorithm identified eight key genes (Figure 6C; Supplementary Table S9), employing a radial basis function kernel with hyperparameters (cost = 0.5, gamma = 0.015625) optimized via 5-fold cross-validation, yielding a minimum cross-validation error rate of 0.28. A Venn diagram showed the key genes identified by both methods: LRRK2, NLRP3, KLRK1, and PLG (Figure 6D). Box plots showed notable disparities in the expression levels of four pivotal neuroinflammation-related genes between the normal and sciatica groups (Figure 6E).

The receiver operating characteristic (ROC) curve showed that a neuroinflammatio-n gene distinguished the normal from the sciatica group (Supplementary Table S10). The area under the curve (AUC) for each individual gene was derived from univariate ROC analysis using its normalized expression value as a continuous predictor. The area under the curve (AUC) and its 95% confidence interval for each gene were as follows: LRRK2: 0.827 (95% CI: 0.704-0.950) (Figure 7A); NLRP3: 0.747 (95% CI: 0.607-0.887) (Figure 7B); KLRK1: 0.800 (95% CI: 0.667-0.933) (Figure 7C); PLG: 0.706 (95% CI: 0.558-0.853) (Figure 7D). Following our two-stage analytical workflow, data from GSE124272 and GSE150408 were integrated for model development. Consensus hub genes (LRRK2, NLRP3, KLRK1, PLG) were identified from this combined dataset through the aforementioned LASSO and SVM-RFE analyses with cross-validation. A final diagnostic model was then constructed using only these four hub genes. The ROC curve of the final model was analyzed together with the aforementioned multiple biomarkers (Supplementary Table S11). This model achieved an AUC of 0.909 (CI = 0.828-0.989) in the development phase (Figure 7E). For the final, unbiased performance assessment, this diagnostic model, constructed using only the identified hub genes, was applied to the pre-specified, hold-out GSE150408 validation cohort. The model achieved an AUC of 0.907 (95% CI = 0.800-1.000) on this internal validation set (Figure 7F). A comprehensive set of performance metrics for this final diagnostic model, including accuracy, sensitivity, specificity, precision, F1-score, and the Matthews Correlation Coefficient (MCC) calculated at the optimal cut-off value, is provided in Supplementary Table S12.

The contributions of the four key neuroinflammation genes to disease risk were evaluated using nomograms (Figure 8A). We also evaluated the clinical utility and prediction accuracy of the models using decision curve analysis and calibration curve analysis. The findings from the decision curve analysis revealed that employing the key genes for risk assessment offers considerable clinical advantages (Figure 8B). Furthermore, the calibration curve demonstrated a strong agreement between the predicted probabilities and the observed outcomes, indicating good model calibration when tested on the validation dataset (Figure 8C; Supplementary Table S13). Using GSEA, associations were identified between groups with high or low expression of four key genes and their corresponding pathological pathways. Among the key genes, LRRK2 is associated with the WP microglia pathogen phagocytosis pathway (NES = 2.201, FDR <0.001) and the KEGG nod-like receptor signaling pathway (NES = 1.969, FDR <0.001) (Figure 9A). NLRP3 is linked to the KEGG nod-like receptor signaling pathway (NES = 1.890, FDR <0.001) as well as the WP toll-like receptor signaling pathway related to MYD88 (NES = 2.047, FDR <0.001) (Figure 9B). KLRK1 shows associations with the REACTOME cellular response to hypoxia (NES = 1.592, FDR <0.001) and REACTOME interferon signaling (NES = 1.867, FDR <0.001) (Figure 9C). Finally, PLG is connected to the NABA core (NES = 1.609, FDR <0.001) and the KEGG neuroactive ligand-receptor interaction pathway (NES = 1.666, FDR <0.001) (Figure 9D).

This study examines changes in the immune cell composition of patients with sciatica using gene expression analysis specific to individual cell types (Figure 10A). We analyzed the immune cell composition between the control and sciatica groups. This analysis revealed numerous statistically meaningful distinctions, particularly in T cells gamma-delta, macrophages M0, and neutrophils. The proportion of M0 macrophages and neutrophils in the sciatica group was significantly greater than in the normal group (p < 0.05). Additionally, the proportion of gamma-delta T cells in the sciatica group decreased significantly compared to the normal group (p < 0.001) (Figure 10B). Spearman correlation analysis revealed a strong positive correlation between activated CD4 memory T cells and various cell types, including activated mast cells, dendritic cells (both activated and resting), naive CD4 T cells, CD8 T cells, plasma cells, and naive B cells (Figure 10C). LRRK2, NLRP3, KLRK1, and PLG correlate significantly with the prevalence of various immune cell types. LRRK2 exhibited a positive correlation with neutrophils and a negative correlation with γδ T cells, CD8 T cells. NLRP3 showed a positive correlation with neutrophils. KLRK1 demonstrated a positive correlation with γδ T cells, resting dendritic cells and macrophages, activated mast cells, and neutrophils. Similarly, PLG was positively correlated with γδ T cells and resting dendritic cells (Figure 10D). The heatmap shows the distribution of immune cells in the normal and sciatica groups (Figure 10E). To address potential limitations in neurological disease applications, we validated CIBERSORTx estimates using MCP-counter analysis. Both methods consistently identified elevated neutrophil infiltration in patients with sciatica, confirming the robustness of our findings across deconvolution platforms (Supplementary Figure S2).

A network of neuroinflammation-related genes involving mRNAs and miRNAs was created, consisting of four mRNAs and sixty miRNAs, with NLRP3 targeted by pro-inflammatory miR-22-3p and miR-223-3p, and KLRK1 showing the most extensive regulatory connectivity, suggesting complex post-transcriptional modulation of the neuroinflammatory signature (Figure 11A). The mRNA-TF network of neuroinflammation-related hub genes, comprising four mRNAs and fifty-nine TFs, was constructed using the ChIPBase v3.0 database, revealing distinct regulatory clusters wherein STAT3 and GATA2 targeted KLRK1, SPI1 and IRF1 regulated NLRP3, and HNF4A exclusively bound PLG, implicating coordinated activation of myeloid inflammation and interferon signaling pathways (Figure 11B). The DGIdb version 3.0.2 (https://www.dgidb.org) predicted potential drugs or small molecule compounds associated with four mRNAs and thirty-two small molecule drugs, notably identifying the FDA-approved IL-1 receptor antagonist Anakinra targeting NLRP3, multi-kinase inhibitors (Vandetanib, Entrectinib) targeting LRRK2, and fibrinolytic agents (Alteplase, Streptokinase) converging on PLG, thereby highlighting actionable therapeutic targets for sciatica intervention (Figure 11C).

In addition, we measured the expression levels of LRRK2, NLRP3, KLRK1, and PLG in patients with sciatica and normal individuals using qRT-PCR. The study found that LRRK2 (p < 0.05, Figure 12A) and NLRP3 were significantly upregulated in sciatica samples compared to standard samples (p < 0.001, Figure 12B), while KLRK1 was significantly downregulated (p < 0.001, Figure 12C). The expression of PLG tends to increase, but not significantly (p > 0.05, Figure 12D).

Although this study presents important findings, it also has several notable limitations. First, the initial gene screening used a p-value threshold rather than a stricter FDR correction to retain an adequate number of candidate genes for downstream machine learning analysis. While this approach was necessary for our hypothesis-generating design and for the robustness of the final model, which was confirmed through rigorous machine learning regularization techniques, it may have increased the potential for false positives in the initial selection phase. Second, the study did not provide a thorough analysis of sciatica-related tissues, and without direct samples, our understanding of the genetic mechanisms underlying neuroinflammation is limited. Third, the machine learning validation employed an internal “leave-one-dataset-out” strategy: data from GSE124272 and GSE150408 were integrated for model development and feature selection via cross-validation, after which the entire, pre-specified GSE150408 cohort was held out exclusively for final performance estimation (AUC = 0.907). While this pre-defined, two-stage workflow rigorously separates model tuning from final evaluation to mitigate optimism bias, it remains an internal validation due to the use of the same data source for both purposes. The scarcity of additional public transcriptomic datasets for sciatica precluded fully independent external validation, which represents a limitation for the generalizability of the biomarker signature. Future studies should increase sample sizes and conduct multicenter research. Fourth, due to the heterogeneity of sciatica and the absence of clinical data, it was impossible to evaluate the relationship between risk indicators and patient stratification. Future research should incorporate detailed clinical data for a comprehensive analysis and interpretation. Fifth, the symptoms of patients and the expression of biomarkers vary, and future research needs to explore the characteristics of different subgroups to identify more specific biomarkers. Sixth, while the use of β-actin as a single reference gene in qRT-PCR experiments was substantiated by rigorous validation demonstrating stability across experimental conditions (including non-significant Cq value differences between groups and minimal technical variability), this approach may not represent the optimal strategy compared to employing multiple reference genes. We have addressed this methodological consideration in the discussion and will incorporate additional validated reference genes in future investigations to align with best practices. Seventh, regarding the assessment of batch effects, the Principal Component Analysis (PCA)-based metric (variance proportion on PC1+PC2) used in this study, while transparent and interpretable for our purpose, may not fully decompose variance components as robustly as advanced methods like Principal Variance Component Analysis (PVCA). Future studies employing multi-batch designs would benefit from implementing such integrated approaches for more granular analysis. Lastly, while this study primarily focused on gene expression analysis, future research should validate the functions of key genes using cell-based experiments and animal models. In summary, future research should increase sample sizes and combine histological and clinical data to more comprehensively elucidate the pathological mechanisms of sciatica.