The two data sets, GSE38484 and GSE27383, were obtained from the GEO database ¹ (22) for analysis. The GSE38484 data set, consisting of 106 SCZ patients and 96 normal control whole blood samples, was used as a training set. The GSE27383 data set, containing peripheral blood mononuclear cell (PBMC) samples from 43 SCZ patients and 29 normal HC, was used as a validation set. Furthermore, 16 genes associated with PMRGs were identified in the GeneCards database ² (23) through the keyword “Propionate metabolism”.

Differential analysis was conducted using R software v4.2.2 and the limma software package (v1.38.0) (24). Genes meeting the criteria of |log2FC| > 0.5 and P < 0.05 were identified as differentially expressed genes (DEGs). The ggplot2 software package (v3.4.1) (25) was used to generate volcano plots and heat maps.

The GSE38484 dataset was used to analyze the expression differences of PMRGs between samples with schizophrenia (SCZ) and control samples. Boxplots were generated using the ggplot2 software package (v3.4.1). The difference in PMRG scores between SCZ samples and HC was computed using the GSVA software package (v1.42.0) (26). Significant differences in PMRG scores between SCZ and control samples were assessed using the Wilcoxon test. Spearman correlation analysis was conducted to determine the correlation between DEGs and PMRGs scores, with a correlation coefficient |r| > 0.5 and P < 0.05 as the criteria for identifying DE-PMRGs.

The distribution of DE-PMRGs on the chromosome was depicted using the Circos visualization tool. The relationship between these genes was then analyzed using Spearman correlation and visualized through a heat map. GO and KEGG pathway enrichment analysis was then conducted on DE-PMRGs using the clusterProfiler software package (v4.7.1.3) (27).

Lasso regression analysis was conducted on the candidate genes using the glmnet software package (v4.1.4) (28). In the Lasso model, the cross-validation method is employed to select the optimal regularization parameter, λ. Specifically, by minimizing the mean squared error (MSE), the optimal value of λ was determined in this study to be 0.0127193 (logλ =-1.8955). The Lasso coefficient spectrum and cross-validation error plots were generated to identify the key genes, which are those with regression coefficients not penalized to 0.

Box plots were generated to display the expression of key genes in SCZ samples and control samples from the training set GSE38484 and validation set GSE27383. Genes identified by their expression levels were noted as potential biomarkers.

Binomial logistic regression analysis was used to develop a diagnostic model using biomarkers from the training set GSE38484 and validation set GSE27383. The ROC package pROC was used to generate an ROC curve to evaluate the model’s efficacy in differentiating SCZ samples from control samples.

In the GSE38484 training set, biomarkers are used to predict disease status as the outcome event. The R package “rms” was used to develop a nomogram model based on these biomarkers, and a calibration curve was generated to assess the predictive performance of the nomogram model. The Hosmer-Lemeshow test (HL test) serves as an indicator of model fit by assessing the discrepancy between predicted and actual values. A p-value greater than 0.05 indicates that the model passes the HL test, suggesting no significant difference between the predicted and actual values. Conversely, a p-value less than 0.05 implies that the model fails the HL test, highlighting a significant discrepancy and indicating poor model fit. Furthermore, a Mean Absolute Error (MAE) of more than 0.05 indicates a small error between actual and predicted disease risks, reflecting the high precision of the nomogram model in predicting disease risk. Additionally, DCA curve analysis was conducted using the R package “rmda” to evaluate the clinical utility of the nomogram. To further assess the effectiveness of the nomogram, the ROC curve was plotted using the R package “PRROC” (29) and the model’s AUC was calculated to evaluate the model’s effectiveness.

The GeneMANIA database ³ (30) was used to construct and analyze the interaction network between biomarkers and their co-expressed genes. Correlations between biomarkers and all genes in the GSE38484 dataset were calculated, and the genes were ranked based on the correlation coefficient. GSEA enrichment analysis was subsequently conducted using the clusterProfiler software package (v4.7.1.3) with reference to the c2.cp.kegg.v2023.1.Hs.symbols.gmt gene set and c5.go.v2023.2 in the MSigDB database ⁴ (31). The top three positive and negative pathways with enrichment scores were selected for visualization based on thresholds of |NES| > 1 and P < 0.05.

The GSVA software package (v1.42.0) was used to assess immune cell infiltration in all samples from the training set GSE38484 using the ssGSEA algorithm. A box plot was generated to illustrate the distribution proportion of 28 immune cell types in the samples. The Wilcoxon test was used to compare differences in immune cell infiltration among the samples. A significance level of P < 0.05 was established for screening, and a box plot was created to visualize the immune scores between samples from individuals with SCZ and control samples. Spearman correlation analysis was used to investigate the relationship between biomarkers and differential immune infiltrating cells. A correlation with |r| > 0.4 and P < 0.05 was deemed statistically significant.

The TarBase database ⁵ (32) was used to predict the miRNAs associated with the biomarkers, while the Starbase database ⁶ (33) was used to predict the lncRNAs corresponding to the miRNAs. Subsequently, a ceRNA network consisting of lncRNA-miRNA-mRNA was established. Moreover, transcription factors (TFs) for the biomarkers were predicted using the ChEA3 database ⁷ (34), where TFs with scores below 500 were selected, creating a TF-mRNA network. The visualization of both the ceRNA network and TF-mRNA network was done using Cytoscape software.

The DrugBank database ⁸ (35) was used to identify potential drugs targeting biomarkers, while Cytoscape software was used to visualize the mRNA-drug network. The 3D structure of the drug was obtained from the PubChem database ⁹ (36), and the protein structure of the biomarker was retrieved from the PDB database ¹⁰ (37). Molecular docking was then performed.

R version 4.2.2 was utilized for all statistical tests. Figure panels were pieced together by EdrawMax. The significance of the correlation between the two groups was assessed through Spearman correlation analysis. Among patients and healthy volunteers, demographic and clinical behavioral data were compared using t-tests or chi-square tests, with a significance threshold set at P < 0.05. Data visualization was carried out using GraphPad Prism v.9.5.1 in conjunction with R version 4.2.2.

The differential expression analysis results of the GSE38484 dataset revealed 27 differential genes between SCZ and normal HC, with 26 up-regulated genes and 1 down-regulated gene. The volcano plot ( Figure 1A ) displays the 10 most significant up-regulated and down-regulated genes, including RPS15A, RPL9, HINT1, EVI2A, RPS3A, COX7C, COMMD6, RPL17, EIF1AY, and RPS4Y1 for up-regulated genes, and NRGN for the down-regulated gene. Furthermore, the heatmap ( Figure 1B ) displays the expression levels of the 27 DEGs.

The boxplot ( Figure 1C ) illustrates the expression variances of PMRGs between samples from individuals with SCZ and HC samples. Among these, 5 PMRGs exhibited significantly higher levels in the control samples compared to the SCZ samples, specifically ACLY (P < 0.01), FBP1 (P < 0.0001), HCFC1 (P < 0.0001), MMAB (P < 0.01), and PCCA (P < 0.001).

The Wilcoxon test evaluated the GSVA scores of differentially expressed PMRGs between the disease group and the HC group within the GSE38484 dataset, revealing a significantly lower score in SCZ compared to the normal HC group (P < 0.05) ( Figure 1D ).

Spearman correlation analysis was conducted to assess the relationship between the DEGs and PMRG scores. Following the screening criteria of correlation coefficient |r| > 0.5, P < 0.05, a total of 11 DE-PMRGs were identified, including ARGLU1, CAPZA2, EVI2A, LY96, NRGN, RPL17, RPL9, RPLP0, S100A8, SUZ12, and TMEM123 ( Figure 1E ).

The Circos heat map illustrates the distribution of 11 DE-PMRGs across different chromosomes. Specifically, S100A8 is situated on chromosome 1, RPL9 on chromosome 4, CAPZA2 on chromosome 7, LY96 on chromosome 8, TMEM123 and NRGN on chromosome 11, RPLP0 on chromosome 12, ARGLU1 on chromosome 13, EVI2A and SUZ12 on chromosome 17, and RPL17 on chromosome 18 ( Figure 2A ).

Spearman correlation analysis was conducted to assess the correlation between DE-PMRGs. The heat map ( Figure 2B ) revealed a significant correlation among all genes (P < 0.05). Notably, NRGN exhibited a significant negative correlation with other candidate genes (r <-0.4, P < 0.05), while the rest of DE-PMRGs showed a significant positive correlation (r > 0.4, P < 0.05).

GO and KEGG enrichment analyses were conducted on DE-PMRGs, resulting in a total of 109 outcomes from GO enrichment. The screening criteria used was P < 0.05, with 69 BP (biological processes), 21 CC (cell components), and 19 MF (molecular functions) being enriched. The enriched results were arranged in ascending order based on the P value, and the top 10 gene functions were presented ( Figure 2C ). Regarding biological processes, DE-PMRGs exhibited significant associations with cytoplasmic translation, positive regulation of lipopolysaccharide-mediated signaling pathways, response to glycoproteins, glial cell differentiation, and positive regulation of cell size. In terms of cell components, DE-PMRGs were notably linked to cytoplasmic large ribosomal subunits, cytoplasmic ribosomes, large ribosomal subunits, ribosomal subunits, and ribosomes. DE-PMRGs were found to bind to TOLL-like receptors, ribosomal structural components, ribosomal RNA, RAGE receptors, and lncRNAs for molecular functions. Furthermore, KEGG functional enrichment analysis identified 13 pathways, with the top 10 enriched pathways including ribosome, novel coronavirus, pertussis, polycomb protein inhibitory complex, IL-17 signaling pathway, NF-kappa B signaling pathway, TOLL-like receptor signaling pathway, toxoplasmosis, alcoholic liver disease, and motor protein pathway. The most significant enrichment results for pathways are illustrated in Figure 2D .

Lasso regression analysis was conducted on 11 DE-PMRGs ( Figures 3A, B ). 5 specific genes, namely LY96, NRGN, RPLP0, S100A8, and TMEM123, were identified as key genes as their regression coefficients were not penalized to 0.

The box plot shows the expression levels of key genes in SCZ samples and control samples in the training set GSE38484 ( Figure 3C ) and the validation set GSE27383 ( Figure 3D ). LY96 and TMEM123 were expressed similarly and significantly differently in these two data sets and were, therefore, identified as biomarkers.

A biomarker-based diagnostic model was constructed in the two data sets using binomial logistic regression analysis, and a ROC curve plot was drawn to evaluate and verify the effectiveness of the model in distinguishing SCZ and control samples. The results show that the AUC of the diagnostic model in the training set GSE38484 is 0.735 ( Figure 4A ), and the AUC in the validation set GSE27383 is 0.847 ( Figure 4B ). In both data sets, the AUC values of the logistic regression models constructed by LY96 and TMEM123 both exceeded 0.7, indicating a good diagnostic effect.

A nomogram ( Figure 4C ) was constructed to predict the probability of SCZ for samples in the training set GSE38484. The Hosmer-Lemeshow test was conducted to determine the dispersion between the predicted and true values. As shown in Figure 4D , the P value of the HL test is 0.568 (> 0.05), indicating no conspicuous difference between the predicted value and the true value. Moreover, the mean absolute error (MAE) is 0.02 (< 0.05), indicating that the error between the actual disease risk and the predicted disease risk is small and that the nomogram model has high accuracy in predicting the disease risk of the sample. In the DCA curve ( Figure 4E ), the nomogram model’s benefit rate, represented by the red line, significantly outperforms both the diagonal line (All) and the horizontal line (None). This indicates that, in comparison to the strategies of either treating all patients or none, the nomogram model offers additional benefits and demonstrates superior performance. Furthermore, in the ROC curve ( Figure 4F ), the model’s AUC exceeds 0.7, suggesting that the nomogram possesses a considerable level of predictive accuracy.

The GeneMANIA database was queried, and 20 genes interacting with the biomarkers LY96 and TMEM123 were obtained. Figure 5A shows the GeneMANIA network of biomarkers, including a total of 7 interactions between biomarkers and genes. Physical interactions account for the highest proportion (77.64%), followed by co-expression (8.01%). Moreover, among the biomarkers and their interacting genes, the top five related functions are the toll-like receptor signaling pathway, response to molecule of bacterial origin, pattern recognition receptor signaling pathway, regulation of pattern recognition receptor signaling pathway, and regulation of toll-like receptor signaling pathway.

The GSEA enrichment analysis results of LY96 and TMEM123 showed that KEGG of LY96 was enriched in 30 pathways and GO was enriched in 46 pathways; KEGG of TMEM123 was enriched in 18 pathways, and GO was enriched in 34 pathways. According to the threshold of |NES| > 1 and P < 0.05, the top three pathways with positive and negative enrichment scores were selected for picture drawing. KEGG enrichment analysis of LY96 found pathways such as ribosomes, autophagy regulation, Leishmania infection, phenylalanine metabolism, tricarboxylic acid cycle, and circadian rhythm disorders ( Figure 5B ). The GO enrichment analysis of LY96 includes pathways such as response to nitrosative stress, cytoplasmic small ribosome subunits, cytoplasmic ribosomes, regulation of nucleotide excision repair, npBAF complex, bBAF complex, etc ( Figure 5C ). KEGG enrichment analysis of TMEM123 revealed pathways such as protein export, ribosomes, basal transcription factors, base excision repair, glycerophospholipid metabolism, porphyrin, and chlorophyll metabolism ( Figure 5D ). GO enrichment analysis of TMEM123 involves pathways such as mitotic sister chromatid cohesion, cytoplasmic small ribosome subunits, GDP binding, actin filament binding, semaphorin plexin signaling pathway, and structural components of muscle ( Figure 5E ).

Figure 6A illustrates the infiltration of immune cells in all samples from the training set GSE38484. The immune cell type with the highest relative abundance is myeloid-derived suppressor cells (MDSC), while the immune cell type with the lowest relative abundance is Th2 cells. Figure 6B presents the differences in immune cell infiltration between SCZ and control samples. Significant variations were observed in 12 types of immune infiltration cells between the disease group and the HC group, including activated CD4 T cells, CD56 bright natural killer cells, CD56 dark natural killer cells, central memory CD8 T cells, effector memory CD4 T cells, effector memory CD8 T cells, eosinophils, γδ T cells, immature dendritic cells, natural killer T cells, follicular helper T cells, and type 1 T helper cells. Then, Spearman correlation analysis was used to calculate the correlation between the biomarkers LY96, TMEM123, and differential immune infiltrating cells (correlation |r| > 0.4, P < 0.05 results were considered as significant correlation results).As depicted in Figures 6C, D , Type 1 T helper cell, Activated CD4 T cell, LY96, and TMEM123 all showed a significant positive correlation(r > 0.4, P < 0.05); LY96, CD56dim natural killer cell, and Effector memory CD8 T cell showed a significant positive correlation. Negative correlation(r <-0.4, P < 0.05), TMEM123 and Central memory CD8 T cell, T follicular helper cell, and Natural killer T cell showed significant negative correlation(r <-0.4, P < 0.05).

The miRNAs of 72 biomarkers were predicted using TarBase. Among them, LY96 predicted 14 miRNAs, and TMEM123 predicted 67 miRNAs. The lncRNAs corresponding to miRNAs were predicted through Starbase. Among them, 72 miRNAs predicted a total of 202 relationships with corresponding lncRNAs and a total of 34 non-repeating lncRNAs. The construction of the lncRNA-miRNA-mRNA ceRNA network ( Figure 7A ) aids in understanding the gene expression regulatory mechanisms of biomarkers.

TFs for biomarkers were predicted using the ChEA3 database. LY96 predicted 42 TFs, and TMEM123 predicted 73 TFs, including 11 common TF transcription factors. A TF-biomarker regulatory network was constructed ( Figure 7B ).

Potential drugs targeting biomarkers LY96 and TMEM123 were searched for through the DrugBank database. LY96 predicted four drugs: Lauric acid, (R)-3-hydroxytetradecanoic acid, Morphine, and Myristic acid ( Figure 8A ). The 3D structures of these drugs and the protein structure of LY96 (PDB ID: 2E56) were obtained in PubChem and PDB databases, and molecular docking was performed on LY96 and the four drugs, respectively ( Figure 8B , Supplementary Figure 2 ). A docking score below-5 kcal/mol indicates that the selected drug has a high binding affinity to the target. In Table 3 , the docking scores of the four drugs and LY96 are all lower than-5 kcal/mol, indicating that small molecule compounds and key gene proteins have high binding affinity. Among them, the one with the lowest docking score was Morphine. Figure 8B shows the molecular docking of LY96 and Morphine.

The expression levels of the LY96 and TMEM123 genes in whole blood samples from patients with SCZ were assessed using quantitative polymerase chain reaction (qPCR). As anticipated, the results indicated that the expression levels of LY96 and TMEM123 in the whole blood of SCZ patients were significantly elevated compared to those in HC ( Figure 9 ).