As depicted in Figure 1, the comprehensive workflow adopted in this study is outlined. We downloaded epilepsy data from the gene expression omnibus (GEO) (Barrett et al., 2012) (Supplementary Table 1). GSE168375 includes 31 epilepsy patient brain tissue samples and 12 normal brain tissue samples as controls (Krawczyk et al., 2022). GSE186334 includes 24 epilepsy patient cortical samples and 12 healthy control cortical samples (Gomes-Duarte et al., 2022). GSE140393 validation set consists of 12 epilepsy patient cortical samples and 9 healthy control cortical samples (Pai et al., 2022). Additionally, we downloaded a set of 580 genes associated with apoptosis from PMID36341760 (Zou et al., 2022) (Supplementary Table 2).

GSE168375 Dataset: Control Group: 12 samples, average age 35.25 (3 males, 9 females). Epilepsy Group: 31 samples (1 missing gender/age), average age 6.7 (13 males, 17 females).

GSE186334 Dataset: Epilepsy Group: 10 samples with available clinical information. Average age for patients E5, E6, E7: 56.667 (2 females, 1 male).

GSE140393 Dataset: Control Group: 9 samples from 3 patients, average age 35.3 (2 males, 1 female). Epilepsy Group: 21 samples from 5 patients, average age 24.8 (4 males, 1 female).

We used the “sva” package (v3.40.0) in R to perform batch correction on the GSE168375 and GSE186334 datasets (Leek et al., 2012). Subsequently, principal component analysis (PCA) was conducted. Differential analysis was performed on the integrated epilepsy dataset using the “limma” package (v3.48.3) to identify DEGs between the epilepsy and control groups (Ritchie et al., 2015). Genes with logFC (fold change, FC) > 0.58 and p-value < 0.05 were considered upregulated DEGs, while genes with logFC < −0.58 and p-value < 0.05 were considered downregulated DEGs. The correlation analysis of DEARGs was performed using the “psych” package (v2.3.6) in R to calculate the Spearman correlation between gene expressions. Gene positions on chromosomes were obtained using the “biomaRt” package (v2.54.1) and visualized using the genes package (v0.5.0) (Kasprzyk, 2011).

The “clusterProfiler” package (v4.7.1.3) is used for performing gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses (Yu et al., 2012). The enrichment results are visualized using the “GOplot” package (v1.0.2) and “enrichplot” package (v1.18.4) (Walter et al., 2015; Wu et al., 2021). The “pathview” package (v1.40.0) in R is used for visualizing the KEGG pathways (Luo and Brouwer, 2013). Gene set enrichment analysis (GSEA) is a computational method used to determine whether a predefined set of genes shows statistically significant differences between two biological states (Subramanian et al., 2005). We downloaded the reference gene set “c2.cp.kegg.v7.4.entrez.gmt” from the MSigDB database (Liberzon et al., 2015). The “clusterProfiler” package (v4.0.5) includes GSEA methods for conducting enrichment analysis and visualization of the dataset (Wu et al., 2021). Gene set variation analysis (GSVA) is used to analyze the variation of gene sets (Subramanian et al., 2005). We obtained the reference gene set “h.all.v7.4.symbols.gmt” from the MSigDB database to perform GSVA analysis on different groups of the integrated GEO dataset. The “GSVA” package (Hanzelmann et al., 2013) (v1.40.0) converts the expression matrix into a pathway enrichment score matrix and uses the “lmFit” function from the “limma” package (v3.48.0) to identify different pathways and calculate p-values (Ritchie et al., 2015).

We used the Lasso regression method with the glmnet R package (v4.1-2) to screen all DEARGs (Friedman et al., 2010). Then, we employed the SVM-RFE algorithm with the caret R package (v6.0.94) to further select features associated with epilepsy-related apoptosis genes. The intersection of the features selected by Lasso and SVM-RFE was considered as the set of genes for constructing the diagnostic model. We performed a Friends analysis using the R package GOSemSim (v2.24.0) (Yu, 2020). Using logistic regression based on the selected feature genes, we constructed the diagnostic model. The diagnostic score was calculated using the gene expression levels and the coefficients obtained from the multiple regression model as follows:

Furthermore, we divided the combined epilepsy dataset into a high-score diagnostic group and a low-score diagnostic group. We used the “pROC” R package (v1.18.0) to plot the ROC curve for the diagnostic model and epilepsy status (Robin et al., 2011). The calibration plot was generated using the rms R package (v6.2-0) to assess the accuracy and discriminative ability of the diagnostic model. Finally, the decision curve analysis (DCA) plot was created using the “ggDCA” R package (v1.1) to evaluate the accuracy and discriminative ability of the logistic regression model (Vickers and Elkin, 2006).

We uploaded the expression profile data of the epilepsy dataset to the CIBERSORTx website and used the provided LM22 gene set to calculate the abundance of 22 immune cell types in each patient (Wu et al., 2021). The correlation between genes related to epilepsy in different groups of immune cells was calculated using the Spearman algorithm, and the correlation heatmap was generated using the R package ggplot2 (v3.3.6).

We conducted further predictive research on the diagnostic model genes using the miRNet database (Chang et al., 2020). We downloaded information on small-molecule drug targets from DrugBank and used this information to predict small-molecule drugs that could interact with the diagnostic model genes (Wishart et al., 2008). We utilized Cytoscape (v3.8.2) for visualization and network analysis (Shannon et al., 2003).

We downloaded the 2D or 3D structures of small molecule drugs from the PubChem database and used Chem3D software to convert the 2D structures into 3D structures (Kim et al., 2021; Bateman et al., 2023). Next, we searched the UniProt and PDB databases for human receptor protein structures corresponding to the genes of interest (Burley et al., 2017). With PyMOL software, we visualized the protein structures (Seeliger and de Groot, 2010). Finally, we used Autodock Vina software to identify the active binding sites involved in the interaction between the ligands and receptors (Seeliger and de Groot, 2010).

All data calculations and statistical analyses were conducted utilizing the R programming language (version 4.2.3). For multiple testing corrections, we utilized the Benjamini–Hochberg method with false discovery rate adjustment. For comparisons between two groups of continuous variables, we estimated the statistical significance of normally distributed variables using an independent Student’s t-test. Differences in non-normally distributed variables were analyzed with the Wilcoxon test. We used Spearman correlation analysis to calculate correlation coefficients between different molecules. All p-values were two-sided, and statistical significance was defined as p < 0.05.

GSE168375 and GSE186334 underwent batch correction, and the pre- and post-batch effect removal data sets were compared using a distribution boxplot and PCA plot (Figure 2).

We intersected the DEGs between the epilepsy and control groups with apoptosis-related genes, resulting in 16 DEARGs, 10 of which were upregulated and 6 were downregulated. Figure 3A displays a volcano plot of the DEGs, Figure 3B depicts a heatmap of the DEGs, and Figure 3C shows boxplots of the DEARGs. The volcano plot, heatmap, and boxplots illustrate the expression patterns of genes, including S100A8 and S100A9, which experienced downregulation in the epilepsy group, whereas PAWR, AR, CD38, and others were upregulated. The Spearman correlation coefficients were calculated amongst the DEARGs (Figure 3D). It was discovered that the expression of HSPA1B and HSPA1A exhibited the highest correlation (r = 0.85, p < 2e-16), while the expression of IL1B and TNF showed positive correlation (r = 0.82, p < 2e-16). Figure 3E illustrates the chromosomal positions of the 16 DEARGs. The figure indicates that AR is situated on the X chromosome whereas S100A8 and S100A9 are located on chromosome 1.

The results (Supplementary Table 3) revealed that in the epilepsy group, the enriched GO pathways were primarily associated with BP, including GO:2001233 (regulation of apoptotic signaling pathway), GO:0051092 (positive regulation of NF-kappaB transcription factor activity), GO:2001237 (negative regulation of extrinsic apoptotic signaling pathway), among others (Figure 4A). In terms of CC, the enriched GO terms included GO:0016235 (aggresome), GO:0016234 (inclusion body), GO:0030667 (secretory granule membrane), and others (Figure 4B). For MF, the enriched GO terms included GO:0050786 (RAGE receptor binding), GO:0035325 (Toll-like receptor binding), GO:0140776 (protein-containing complex destabilizing activity), and others (Figure 4C). Figure 4D illustrated the enriched KEGG pathways, mainly involving the IL-17 signaling pathway, Antigen processing and presentation, and others. Figure 4E presented a heat map of the enriched KEGG pathways, providing a visual representation of the pathway-gene associations (Supplementary Table 4).

Supplementary Figure 1A illustrated the role of genes in the TNF signaling pathway. In this pathway, the inflammatory factors TNF and IL1B activate TNFR1 and TNFR2, respectively, leading to the activation of downstream pathways that can induce cell apoptosis or perform other functions. Supplementary Figure 1B depicted the specific details of the Antigen processing and presentation pathway. TNF activation triggered a series of subsequent reactions, including the involvement of PA28, which assisted in the functioning of CD8 T cells, NK cells, CD4 T cells, and other immune cells.

Due to the subjective nature of selecting threshold values for DEGs, we also performed GSEA based on gene expression fold changes (Supplementary Table 5). Figure 5A represented the TGF BETA SIGNALING PATHWAY, Figure 5B showed CYTOKINE CYTOKINE RECEPTOR INTERACTION, Figure 5C displayed COMPLEMENT AND COAGULATION CASCADES, Figure 5D illustrated ECM RECEPTOR INTERACTION, Figure 5E depicted HEMATOPOIETIC CELL LINEAGE, and Figure 5F represented RIBOSOME. All of these pathways were enriched in the epilepsy group.

We conducted GSVA based on gene expression profiles (Supplementary Table 6). Supplementary Figure 2A displayed stronger functional enrichment of TGF_BETA_SIGNALING, KRAS_SIGNALING_UP, and EPITHELIAL_MESENCHYMAL_TRANSITION in the epilepsy group, while REACTIVE_OXYGEN_SPECIES_PATHWAY, KRAS_SIGNALING_DN, and HEME_METABOLISM showed weaker functional enrichment in the epilepsy group. We examined the Spearman correlation between the gene with the highest absolute fold change in expression between the two groups (S100A9) and functional pathways. Supplementary Figure 2B presented a scatter plot showing the correlation between the GSVA score of TGF_BETA_SIGNALING and S100A9 gene expression, which exhibited a negative correlation (r = −0.36, p = 0.0013). Supplementary Figure 2C demonstrated a positive correlation between the GSVA score of REACTIVE_OXYGEN_SPECIES_PATHWAY and S100A9 gene expression (r = 0.35, p = 0.0016). Finally, Supplementary Figure 2D showed a negative correlation between the GSVA score of EPITHELIAL_MESENCHYMAL_TRANSITION and S100A9 gene expression (r = −0.32, p = 0.0041).

We employed the Lasso algorithm to select twelve feature genes from sixteen DEARGs (Figures 6A, B). Furthermore, we used the SVM-RFE algorithm to select five genes from the sixteen DEARGs (Figure 6C). Combining the results of both algorithms, we identified five feature genes (CD38, FAIM2, IL1B, PAWR, S100A8) as diagnostic biomarkers for epilepsy grouping (Figure 6D). Using logistic regression, a diagnostic model for epilepsy grouping was constructed with these five key genes. The risk score was calculated as: risk score = 0.596134*exp(CD38)–1.312078*exp(FAIM2) + 0.612583*exp(IL1B) + 1.568740*exp(PAWR)–0.687394*exp(S100A8). The constructed diagnostic model demonstrated high accuracy in diagnosing epilepsy, with an area under curve (AUC) value of 0.916 as shown in Figure 6E.

We used GSE140393 as a validation set and calculated the ROC curve accordingly (Figure 6F). The results showed that the diagnostic model had high accuracy in diagnosing epilepsy in the external validation dataset GSE140393, with an AUC of 0.722. Furthermore, individual ROC curves were generated for each gene separately (Figure 6G). The AUC values ranged from 0.66 to 0.82, indicating that the diagnostic capacity of models based on individual genes was lower than that of the overall 5-gene model. Lastly, we conducted GO analysis to assess the functional importance of the 5 diagnostic genes using Friends analysis (Figure 6H). The results showed that the gene FAIM2 played an important functional role and might be a crucial gene in the context of epilepsy.

We analyzed the impact of each gene on the disease and visualized it as a forest plot (Figure 7A). The results showed that IL1B and PAWR were risk factors. On the other hand, FAIM2 and S100A8 were protective factors. The diagnostic capability of the diagnostic model was assessed using Nomogram analysis (Figure 7B). It was found that PAWR made a significant contribution to the diagnosis of epilepsy. The calibration curve (Figure 7C) is primarily used to assess the accuracy of the diagnostic model. The fitted curve showed a high degree of overlap with the reference curve (dashed line), indicating good predictive performance of the diagnostic model. The decision curve (Figure 7D) represented the net benefit against threshold probabilities. This indicates that the clinical effectiveness of the model was better within this range. Upon observation, it can be concluded that the diagnostic model showed good clinical effectiveness.

We divided the samples into a high-risk diagnostic group and a low-risk diagnostic group based on the median risk score from the diagnostic model and calculated the expression differences between the two groups. Then, we performed GSEA analysis on the high-risk and low-risk diagnostic groups based on the fold changes in gene expression.

Figure 8A showed the GSEA analysis results for the TGF BETA SIGNALING PATHWAY, Figure 8B displayed the GSEA analysis results for ECM RECEPTOR INTERACTION, Figure 8C presented the GSEA analysis results for CYTOKINE CYTOKINE RECEPTOR INTERACTION, Figure 8D illustrated the GSEA analysis results for NOD LIKE RECEPTOR SIGNALING PATHWAY, and Figure 8E represented the GSEA analysis results for FOCAL ADHESION. All five pathways were found to be more enriched in the high-risk diagnostic group. Figure 8F showed the GSEA analysis results for OXIDATIVE PHOSPHORYLATION, which was more enriched in the low-risk diagnostic group.

We performed GSVA analysis on the gene expression profiles. Supplementary Figure 3A presented the results of the functional enrichment analysis, which showed that UV_RESPONSE_DN, TGF_BETA_SIGNALING, KRAS_SIGNALING_UP, and EPITHELIAL_MESENCHYMAL_TRANSITION were functionally enriched to a greater extent in the high-risk diagnostic group, while REACTIVE_OXYGEN_SPECIES_PATHWAY, MYC_TARGETS_V2, KRAS_SIGNALING_DN, and GLYCOLYSIS were more functionally enriched in the low-risk diagnostic group.

Subsequently, we examined the Spearman correlation between the risk diagnostic score and functional pathways. Supplementary Figure 3B showed the scatter plot of the correlation between the GSVA score of GLYCOLYSIS and the risk diagnostic score, revealing a negative correlation (r = −0.3, p = 8.1e-03). Supplementary Figure 3C displayed the negative correlation between the GSVA score of MYC_TARGETS_V2 and the risk diagnostic score (r = −0.39, p = 4.5e-04). Supplementary Figure 3D demonstrated the negative correlation between the GSVA score of REACTIVE_OXYGEN_SPECIES_PATHWAY and the risk diagnostic score (r = −0.45, p = 4.2e-05). Supplementary Figure 3E represented the positive correlation between the GSVA score of TGF_BETA_SIGNALING and the risk diagnostic score (r = −0.41, p = 2.1e-04). Supplementary Figure 3F illustrated the negative correlation between the GSVA score of EPITHELIAL_MESENCHYMAL_TRANSITION and the risk diagnostic score (r = 0.35, p = 1.6e-03).

We compared the differences in immune cell infiltration abundance between the high and low-risk diagnostic groups for each of the 22 immune cell types (Figure 9A). The results showed significant differences in the infiltration abundance of T cells CD4 naive and Eosinophils. T cells CD4 naive had higher immune infiltration in the low-risk diagnostic group, while Eosinophils had higher immune infiltration in the high-risk diagnostic group.

In the high-risk group (Figure 9B), FAIM2 and Monocyte showed positive correlations (r = 0.46, p = 2.6e-02), while they exhibited negative correlation with T cells CD4 naive (r = −0.32, p = 4.7e-02). In the low-risk group (Figure 9C), FAIM2 showed a negative correlation with T cells CD4 naive (r = −0.33, p = 3.8e-02), and IL1B showed a negative correlation with T cells CD4 naive (r = −0.36, p = 2.3e-02).

Figure 10A demonstrated that as the expression of IL1B increased, the proportion of Eosinophils also increased. Figure 10B illustrated that with an increase in PWAR expression, the proportion of Eosinophils also increased. Figure 10C showed that as S100A8 expression increased, the proportion of Eosinophils decreased. Figure 10D displayed that as IL1B expression increased, the proportion of T cells CD4 naive decreased.

The gene-miRNAs network consisted of 164 interaction pairs involving 148 miRNAs (Supplementary Figure 4A and Supplementary Table 8). The gene-TFs network consisted of 22 interaction pairs involving 20 TFs (Supplementary Figure 4B and Supplementary Table 7). The gene-drug small molecule network consisted of 15 interaction pairs involving 15 small molecule drugs (Figure 11 and Supplementary Table 9).

According to the gene-small molecule interaction data, it was determined that two out of the five selected feature genes, IL1B and S100A8, interacted with the small molecules. Following the selection criteria for protein structures (presence of at least one ligand, lower resolution preferred), we identified the 5R8Q structure of IL1B that met the requirements. Subsequently, we performed molecular docking analysis between IL1B (5R8Q) and the small molecules Andrographolide, Binimetinib, Dilmapimod, Etiprednol Dicloacetate, and others.

Figure 12A depicted the binding scenario between IL1B and Andrographolide, where hydrogen bonds with residues LYS-103, MET-148, and ASN-108 were observed, and the predicted binding affinity was −7.2. In Figure 12B, the binding interaction between IL1B and Binimetinib was shown, with a hydrogen bond formed with residue LEU-26, and the predicted binding affinity was −7.2. Figure 12C illustrated the binding of IL1B and Dilmapimod, with hydrogen bonds formed with residues LYS-74 and VAL-163, and the predicted binding affinity was −9.0. Figure 12D displayed the binding between IL1B and IL1B_Etiprednol Dicloacetate, forming a hydrogen bond with residue MET-20, and the predicted binding affinity was −6.1. Meanwhile, we also show the spatial structure of the 5R8Q construct of IL1B (Supplementary Figure 5).