This study retrieved five SCZ-related microarray datasets (GSE12654, GSE21935, GSE17612, GSE53987, GSE38481) from the Gene Expression Omnibus (GEO) database, comprising a total of 224 brain tissue samples (112 SCZ and 112 controls) and 37 whole blood sample (22 SCZ and 15 controls). Detailed dataset information is provided in Table 1 . Quantile normalization was performed using the “limma” package to eliminate technical variability, and the ComBat algorithm was applied to correct inter-platform batch effects. After quality control, GSE12654 and GSE21935 were merged as the training cohort, while GSE17612, GSE53987 and GSE38481 served as independent external validation cohorts. The effectiveness of data integration was validated using boxplots and principal component analysis (PCA) generated by ggplot2.

Differential expression analysis was conducted on the training cohort, GSE12654 and GSE21935 using the “limma” package, with screening criteria set at p-value <0.05 and |logFC| > 0.585 (14, 15). Heatmap and volcano plots of DEGs were generated using “pheatmap” and “ggplot2” packages, respectively.

The RRA algorithm integrates gene ranking information across datasets via a probabilistic model to identify consistently significant DEGs across multiple independent datasets. In this study, RRA was applied to rank up- and down-regulated DEGs from all datasets based on logFC. Aggregated ranking scores were used to compute p-value, and genes with p-value <0.05 and |logFC| > 0.585 were selected as DEGs. RRA analysis was implemented using the “RobustRankAggreg” package.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using the “clusterProfiler” and “org.Hs.eg.db” packages. Significantly enriched GO terms and KEGG pathways were defined as those with p-value <0.05.

PPI networks of DEGs were constructed using the STRING database and visualized with Cytoscape software. Feature DEGs were identified using 10 topological network algorithms (MCC, MNC, etc) via the cytoHubba plugin.

Lasso regression and RF were employed for hub genes selection. Lasso regression is a linear regression method used for feature selection and sparse modeling. It incorporates an L1 regularization term into the objective function to select fewer features, thereby reducing the risk of overfitting. In this study, Lasso regression was implemented using the “glmnet” package with the following parameters: family = binomial, type.measure = class, alpha = 1, and nfold = 10. RF, an ensemble learning method based on decision trees, was employed to capture non-linear relationships and assess feature importance. By constructing multiple decision trees and combining their predictions, RF enhances model accuracy and controls overfitting. The RF algorithm was performed using the “randomForest” package with ntree = 500. Feature importance was evaluated using the Gini coefficient, with a threshold of >2 used for selecting features DEGs.

We utilized the “glm” function from the “glmnet” package to construct a logistic regression model. This function is a standard tool for fitting generalized linear models, particularly suitable for logistic regression analysis in binary classification problems. The expression level of feature DEGs was severed as the independent variable, with the disease diagnosis outcome (SCZ = 1, Control = 0) as the dependent variable. The parameter was set as family = binomial (link=‘logit’). The formula of the model: y= β₀+β₁*X₁+ β ₂*X₂+ β ₃*X₃+⋯+β_i*X_i (β was coefficient, X was the expression level of genes). Diagnostic performance was evaluated via receiver operating characteristic (ROC) curve in the training cohort and external validation cohorts. The nomogram for predicting disease risk was constructed using the “rms” package, and the clinical applicability of the model was assessed through a calibration curve.

single-sample gene set enrichment analysis (ssGSEA) quantified infiltration levels of 28 immune cell subtypes (16, 17). Wilcoxon rank-sum test was used to compare immune cell infiltration between SCZ patients and controls. Spearman’s rank correlation analyzed associations between hub genes and immune cells.

The study workflow was illustrated in Figure 1 . After batch effect correction, showed significant improvements in sample clustering ( Figures 2A-D ). In the training cohort, a total of 29 DEGs were identified, including 16 downregulated and 13 upregulated genes ( Figure 2E , Supplementary Table S1 ). Before RRA analysis, differential expression analyses were conducted separately on the GSE12654 and GSE21935 datasets, and the results were shown by volcano plots ( Figures 3A, B ). After integration using the RRA algorithm, a total of 98 DEGs were identified, comprising 41 downregulated and 57 upregulated genes ( Figure 3C , Supplementary Table S2 ). Taking the intersection of these DEGs, we ultimately identified 22 significant DEGs ( Figure 3D ).

Subsequently, GO and KEGG enrichment analyses were performed on the 22 significant DEGs. The biological processes (BP) with high significance were related to nervous system development and function, including glial cell differentiation, astrocyte differentiation; Cellular component (CC) showed significant enrichment of DEGs in secretory granule lumen, collagen−containing extracellular matrix; Molecular function (MF) revealed significant enrichment of DEGs in carboxylic acid binding and organic acid binding ( Figure 4A ). KEGG pathway analysis further indicated that DEGs were enriched in several key signaling pathways involved in biological processes, such as IL-17, TNF and Hippo signaling pathways ( Figure 4B ).

To systematically identify hub genes associated with disease pathogenesis, we first constructed PPI network to visualize the interactions among 22 DEGs. The results showed that 11 of the 22 DEGs had interactions ( Figure 5A ). Using 10 topological network algorithms to rank genes, the intersection of the top 10 genes from each algorithm was taken, with a total of 7 genes (S100A9, CHI3L1, WWTR1, VGLL1, SERPINA3, S100A8, PVALB) identified as feature DEGs ( Figure 5B ). Subsequently, lasso regression and RF were employed to identify feature DEGs. Lasso regression analysis selected 13 feature DEGs (MAG, VGLL1, S100A8, SPRR1A, ZNF345, S100A9, USH1C, NMU, SH2D2A, GPR45, SERPINA3, ARHGEF5, IGFBP2) ( Figures 5C, D ). while RF identified six feature DEGs (SH2D2A, WWTR1, MAG, VGLL1, KLK6, and S100A9) with Gini coefficients >2 ( Figures 5E, F ). Ultimately, through the intersection analysis of these features DEG subsets, two DEGs (S100A9 and VGLL1) were determined as the optimal hub genes for SCZ ( Figure 6A ).

The diagnostic value of the hub genes in SCZ was further evaluated. In the training cohort, S100A9 and VGLL1 were significantly upregulated in SCZ groups ( Figure 6B ). The area under curve (AUC) values of ROC curve for these genes was 0.702 and 0.694, respectively ( Figure 6C ), indicating their diagnostic potential. Subsequently, a logistic diagnostic model was constructed based on the expression levels of S100A9 and VGLL1, with the formula: y=−1.4815 + 1.2469×VGLL1 + 0.4252×S100A9. The model achieved an AUC of 0.806 in the training cohort, demonstrating good discriminatory ability for SCZ ( Figure 6D ). Subsequently, we performed external validation of the predictive model’s discriminative performance using gene expression profiles of independent validation cohorts (GSE17612, GSE53987, and GSE38481). The results demonstrated that the model achieved AUC values of 0.702, 0.666, and 0.739 in the three external cohorts, respectively, further confirming its diagnostic efficacy and generalizability across independent datasets ( Figures 6E-G ). Additionally, we constructed a nomogram to predict the risk of SCZ ( Figure 7A ), and calibration curve analysis showed high consistency between predicted and actual SCZ risks ( Figure 7B ).

To elucidate the potential molecular mechanisms of hub genes in the development of SCZ, we divided the SCZ samples in the training cohort into high- and low-expression groups based on the median expression of the hub genes and performed differential expression and functional enrichment analyses. Based on the median expression of S100A9, a total of 211 DEGs were identified, including 67 downregulated and 144 upregulated genes ( Figure 8A , Supplementary Table S3 ). GO analysis ( Figure 8B ) indicated that S100A9-related DEGs were primarily involved in BP related to immune system regulation, such as regulation of immune effector process, leukocyte cell−cell adhesion, and cell activation involved in immune response. In terms of CC, DEGs were significantly enriched in external side of the plasma membrane and collagen−containing extracellular matrix. Regarding MF, DEGs were significantly enriched in cytokine binding, immune receptor activity, and cell adhesion mediator activity. Additionally, KEGG pathway analysis showed that DEGs were primarily involved in the PI3K-Akt, HIF-1, and TNF signaling pathways ( Figure 8C ).

Based on the median expression of VGLL1, a total of 54 DEGs were identified, including 30 downregulated and 24 upregulated genes ( Figure 8D , Supplementary Table S4 ). GO analysis showed that the main BP enriched by DEGs were also related to the immune system, such as leukocyte cell-cell adhesion and leukocyte adhesion to vascular endothelial cells; in addition, collagen-containing extracellular matrix, transport vesicles, growth factor activity, and heparin binding were significantly enriched in CC and MF ( Figure 8E ). KEGG pathway analysis further revealed that DEGs were significantly enriched in PI3K-Akt, TNF and IL-17 signaling pathways ( Figure 8F ). These results, similar to those of S100A9-related DEGs, suggest the important regulatory roles of the immune system and signaling pathways in the pathogenesis of SCZ.

Given the enrichment analysis results indicating the involvement of hub genes in immune system responses, we further employed ssGSEA to analyze the infiltration levels of 28 immune cell types in the training cohort. The analysis revealed significant increases in the infiltration levels of five immune cell types in the SCZ group, including CD56bright natural killer cells, MDSCs, mast cells, natural killer cells, and plasmacytoid dendritic cells ( Figure 9A ). Correlation analysis further showed that S100A9 was significantly positively correlated with these immune cells, while VGLL1 was negatively correlated with plasmacytoid dendritic cells ( Figures 9B, C ). Collectively, our results suggest that hub genes may be involved in the pathogenesis of SCZ through the regulation of immune cell infiltration.