From the GEO database, ¹ the schizophrenia data set GSE92538 was selected as the training group and GSE21935 as the test group (13). CellAge database ² was used to retrieve 279 cellular senescence-related genes (14). Figure 1 shows the process flow and the complete information dataset is given in Table 1.

Gene expression profile was used (15) to calculate the MAD (Medium Absolute Development) of all genes in the GSE92538 dataset, and the first 25% of the genes with the least MAD were removed. To eliminate the genes and sample outliers, the GoodSamplesGenes method of the WGCNA package of R software was utilized and further construction of a scale-free co-expression network was also done by WGCNA (the sensitivity was kept 3). Furthermore, modules with a distance of less than 0.25 were combined, and as a result, 13 co-expression modules were obtained. It was important to note that the grey module was regarded as a gene set that was not compatible with any module.

Gene functional enrichment analysis of the above-processed module genes was performed to understand the main expression forms of schizophrenia genes. For gene set functional enrichment analysis, KEGG rest API () was used to acquire gene annotation of the KEGG pathway. Additionally, GO annotations for genes were used from “R” (version 3.1.0) (org.Hs.eg.db) (16). To acquire the findings of the gene set enrichment, they were mapped onto the background set through ClusterProfiler (version 3.14.3) of “R” (17). The GSEA software (version 3.0) was downloaded for GSEA from its website ³ (18). To assess the associated pathways and molecular mechanisms, the samples were separated into two groups and the subsets were downloaded from the Molecular Signatures Database (c2.cp.kegg.v7.4.symbols.gmt) (19). A minimum gene set of 5 and a maximum gene set of 5,000 was determined on the basis of gene expression profiles and phenotypic grouping. p value <0.05 and FDR of <0.1 were considered statistically significant.

The Limma package (linear models for microarray data) (version 3.40.6) (20) of the R software was used for screening DEGs between comparison groups and control groups. In this study, | log2 multiple change (FC) | > 0.5 and p value <0.05 were selected as the criteria for the Limma package to identify differential gene expression (DGE). The related genes were cross-screened to schizophrenia DEG and cell senescence through Venn diagrams (21), and heat maps and volcano plots of cell senescence-related DEG were visualized through sangerBox (22).

In the current study, the gene expression data, its survival time, and survival status were integrated using glmnet (23) and RandomForest (24) of the R software, and regression analysis were carried out using the Random Forest and lasso-cox methods. In addition, to optimize the model, 10-fold cross-validation was set. The neuralnet (25) of “R” was used to construct an artificial neural network for the feature genes obtained by the above methods, to build a high-precision diagnostic model. Further, pROC (26) of “R” was used to conduct ROC analysis to obtain AUC, and the CI tool of the pROC was utilized for the evaluation of AUC and CI (confidence interval) to obtain the final results. The sangerBox software was used for visualization and observation of the expressed characteristic genes in the training set (GSE92538) and test group (GSE21935).

The corrplot of “R” (27) was used to create the heat maps related to invasive immune cells, and Cibersort (28) was utilized to examine the infiltration of immune cells.

The DrugBank database, a free access and comprehensive online database with information about drugs and their targets were used to retrieve the drug-targeting central genes (29). The drugs related to the characteristic genes were recorded and visualized for further correlation analysis.

A cluster tree diagram of schizophrenia and the control group was created using the soft threshold β = 6 (Figure 2A). As a result, 13 GCM (gene co-expression modules) were constructed (Figures 2B,C), and the relationship between GCM and schizophrenia was demonstrated (Figure 2D). Cyan module (124 genes) was found to have the highest correlation with schizophrenia (correlation coefficient = 0.21, p = 1.0e-3). The correlation between module members and gene significance in the cyan module of schizophrenia was also calculated, and a significant positive correlation was observed (r = 0.28) (Figure 2E).

The processed schizophrenia module genes were subjected to functional enrichment analysis. KEGG analysis showed that CSs were mostly present in the “Epstein–Barr virus infection,” “Phagosome” and “Tuberculosis” pathways (Figure 3A). The GO analysis revealed that CGs were primarily found in “vesicle,” “extracellular region,” and “cyclomatic vesicle” sites of cell components (CC; Figure 3B). The key biological processes (BP) of CGs included the immune system process, response to stress and response to an organic substance (Figure 3C). Based on the molecular function (MF) it was revealed that the key components in CGs were “signaling receiver binding” and “identity protein binding” (Figure 3D).

Using the Limma package, 1796 DEGs were identified in the schizophrenia data set (GSE92538), of which 898 were up-regulated and 898 were down-regulated. Schizophrenia DEGs and cell senescence-related genes were cross-screened through Venn diagrams to obtain 23 cell senescence DEGs (Figure 4A). In addition, relevant heat maps and volcano plots were drawn, including the 15 up-regulated and 8 down-regulated genes (Figures 4B,C).

The functional enrichment of genes related to schizophrenia and cell senescence was analyzed. KEGG analysis revealed that CSs were most abundant in the “Epstein–Barr virus infection,” “Toll-like receptor signaling” and “Hepatitis B” pathways (Figure 5A). GO analysis revealed that in terms of CC, CGs were mainly located in “cyclosol,” “nucleoplasma” and “nuclear lumen” (Figure 5B). “Tissue formation,” “positive regulation of gene expression,” and “adaptive process” were among the major biological processes (BP) of CGs (Figure 5C). Based on molecular function (MF), the most significant components of CGs were “protein kinase activity,” “kinase activity,” and “double translated DNA binding” (Figure 5D). By comparing the functional enrichment analysis of schizophrenia, it was revealed that both of them mainly occurred in the pathway related to virus and immunity, which also verified the correlation between them.

LASSO regression analysis was used for the screening of candidate genes and 13 potential biomarkers were identified from these results (Figures 6A,B). Followed by the RF regression analysis for the screening of candidate genes and 13 potential candidate biomarkers were displayed from it (Figure 6C). Then the results from the two aforementioned techniques were cross-analyzed, and finally, six candidate genes (SFN, KDM5B, MYLK, IRF3, IRF7, and ID1) were obtained (Figure 6D). A neutral network was constructed using these six genes, and the results of these confirmed that the six candidate genes could well differentiate the schizophrenia samples from the control samples (Figure 6E). The expression profile analysis of these six genes was evaluated, and results revealed that the candidate genes were considerably different between the control group and the schizophrenia group (Figure 6F).

A nomogram of six potential candidate genes was established (Figure 7A), the CI (95%) and AUC of each item were analyzed, and ROC curves were plotted to determine the specificity and sensitivity of each gene. Nomogram (AUC 0.80, CI 0.87–0.73), SFN (AUC 0.71, CI 0.79–0.63), KDM5B (AUC 0.70, CI 0.78–0.62), MYLK (AUC 0.64, CI 0.73–0.56), IRF3 (AUC 0.65, CI 0.73–0.56), IRF7 (AUC 0.67, CI 0.75–0.59), ID1 (AUC 0.64, CI 0.73–0.55; Figures 7B–H). These findings demonstrated the high diagnostic value of each candidate gene for metabolic syndrome and bipolar illness. In order to further verify the model, a nomograph of the candidate genes in the test group (GSE21935) and ROC curves (AUC 0.83, CI 0.97–0.69) was constructed (Figures 8A,B). The results of the test group showed that the model had certain significance for the diagnosis of schizophrenia. Moreover, the gene selection enrichment analysis (GSEA) on six candidate genes showed its correlation with metabolism (Linoleic acid metabolism, Tryptophan metabolism, Phenylalanine metabolism) (Figures 8C–H).

In the current study, the CIBERSORT algorithm was used to assess the proportion of 22 immune cells in schizophrenia samples and control samples in the training group (GSE92,538; Figure 9A). Followed by the comparison of the immunocyte infiltration of schizophrenic and control samples in the box plot (Figure 9B). The results revealed that the levels of plasma cells, memory B cells, and M2 macrophages in schizophrenic patients were higher, but the levels of regulatory T cells, CD8+ T cells, and follicular helper T cells were lower.

DrugBank database was used to collect data on two drugs (Fostamatinib and Ritodine) that target six candidate genes (Figure 10) and they were the MYLK inhibitors. Fostamatinib was used for the treatment of immune thrombocytopenic purpura (ITP). While Ritodine was an adrenergic β receptor agonist used for the treatment of premature birth.