Gene datasets related to bipolar disorder were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geoprofiles/). The inclusion criteria are set as: expression profiles come from the same sample source and should contain enough sample size that at least fifty samples to ensure accuracy. In addition, the test specimens included should be from humans. Finally, two microarray datasets [GSE12649 and GSE5392] derived from prefrontal cortex were downloaded from GEO (Affymetrix GPL96 platform, Affymetrix Human Genome U133A Array). The GSE12649 comprised 33 bipolar disorder samples and 34 control groups, while the GSE5392 contained 30 bipolar disorder samples and 52 control groups. GSE12649 was selected as the training group and GSE5392 as the test group. The expression matrix file of GSE12649 was normalized through the Limma package (29).

Since no microarray data sets detecting gene expression in brain especially in prefrontal cortex among arteriosclerosis, DisGeNET (https://www.disgenet.org/home/) was used to retrieve 2006 arteriosclerosis-related genes (ASRGs) ( Supplementary Table S1 ). The DisGeNET database is a comprehensive gene-associated information platform that integrates experimentally validated data with information from authoritative repositories and scientific literature to present genes and variants linked to human diseases (30). Given that the GEO database lacks datasets directly related to AS, and considering that atherosclerosis is a consequence of arteriosclerosis (10), we have chosen the atherosclerosis dataset GSE100927 as our test group. The whole process flow is shown in Figure 1 and the complete information dataset is given in Table 1 .

After preparing the data, differential expression analysis was performed using the limma package in R (version 4.3.1) on the GSE12649 dataset. Differentially expressed genes (DEGs) were identified which were differentially expressed between the BD and control groups. The DEG threshold was set as P value <0.05 and |log2FC (fold change) |>0.2. Next, the difference analysis results were presented using the heatmap and volcano plot. In both plots, blue indicated low expression, and red indicated high. Then, the online Venn diagram tool (http://bioinformatics.psb.ugent.be/webtools/Venn/) was used to obtain the differentially expressed arteriosclerosis-related genes (DE-ASRGs) between DEGs and ASRG.

To further understand the function of the DE-ASRGs, we performed Gene Ontology (GO) enrichment analysis (31) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis (32) using the DAVID (https://david.ncifcrf.gov/). GO is a structured way to represent biological functions in terms of core entities and annotate protein biomarkers in the biological process (BP), cellular composition (CC) and molecular function (MF) levels. KEGG enables the correlation of gene catalogues to system functions at the cellular, species and ecosystem levels, facilitating researchers to understand the signaling pathways in which genes are involved. Based on the analysis, we selected the enrichment analysis results with a count greater than 5 and arranged them in descending order of count to identify the top 10.

The DE-ASRGs were analyzed by a search tool for the retrieval of interacting genes/proteins (STRING) (https://www.string-db.org/) to predict the PPI network and to determine the possible relationships between them (confidence level 0.4). The STRING database and Cytoscape software (v3.10.1) completed constructing the PPI network for DE-ASRGs. Gene modules were analyzed and network characteristics of genes were ranked by protein score via the MCODE and CytoHubba plugin, respectively. The CytoHubba plugin was used to score each node gene by 9 randomly selected algorithms, including MNC (Maximum Neighbourhood Component), Degree, MCC (Maximal Clique Centrality), EPC (Edge Percolated Component), Closeness, BottleNeck, Betweenness, Radiality, and Stress. The top 10 genes from each algorithm were used to screen candidate hub genes through the “UpSetR” package.

The WGCNA was utilized for constructing unsigned co-expression networks for the identification of co-expression modules across samples using the “WGCNA” R package (33). WGCNA is a highly efficient and accurate method for analyzing microarray data, and using this approach can help identify genes associated with diseases. Initially, normalized mRNA expression data were used to perform WGCNA to identify gene co-expression and the correlation between gene modules and clinical characteristics (BD compared to control groups). Samples were examined for missing values and then clustered using the average linkage hierarchical method. Afterward, the optimal values of the weighted parameters of the adjacent functions were obtained using the pickSoftThreshold function and were used as soft thresholds for subsequent network construction. Moreover, a topological overlap matrix (TOM) was devised on the basis of an adjacency matrix, and a dynamic tree-cutting algorithm was used for detecting gene modules with a minimum gene group size of 30 and power = 6. The correlation of each module with the BD was calculated, and the module with P<0.05 was defined as the key module. Finally, we took the intersection between DE-ASRGs and WGCNA-derived key module genes to obtain potential hub gene of AS and BD.

To further identify hub genes from the potential candidates, machine learning algorithms were employed. To ensure the repeatability of these algorithms, we set the seed at 2023. The Least Absolute Shrinkage and Selection Operator (LASSO) algorithm was executed using the glmnet package (34), with a tenfold cross-validation was performed to adjust the optimal penalty parameter. The response type was set to binomial, nlambda was set to 100, and alpha was fixed at 1. Moreover, we chose the best lambda value by “lambda.min”. Subsequently, the Random Forest (RF) algorithm (35) was implemented using the randomForest package. We explored the optimal number of random forest trees using cross-validation errors and settled on 500 trees for analysis. The significance of genes in the RF model was evaluated based on mean decrease accuracy and mean decrease Gini. The intersection genes from the top 3 mean decrease accuracy and top 3 mean decrease Gini were identified. Additionally, the support vector machine recursive feature elimination (SVM-RFE) method (36), acting as a vigilant machine learning approach, was employed to select the most relevant variables by eliminating SVM-produced eigenvectors. For the categorization analyses of the screened markers in the BD, the caret package was utilized to conduct the SVM-RFE. The results of SVM-RFE were visualized, and through ten-fold cross-validation, the blue point represented the maximum classification accuracy, highlighting the corresponding gene sets as the most effective diagnostic markers. Only the names of valuable genes used in each machine learning model were included, irrespective of the specific learning model. The hub genes were determined as the intersection genes of valuable genes obtained from the three machine learning algorithms.

SVMs have been shown to exhibit robust classification performance and high accuracy when utilizing feature selection genes in the analysis of microarray expression profiles (37). In this study, SVM was used to evaluate the predictive capability of hub genes as feature genes as follows: 1) the BD datasets including GSE12649 and GSE5392 were separately divided into training (70%) and test (30%) sets; 2) SVM modeling was conducted using gene expression data of hub genes from the training set; 3) optimal SVM hyperparameters were determined through ten-fold cross-validation; 4) model accuracy was assessed using the test set.

To validate the expression of two hub genes in BD and AS, the t-test was employed to analyze their expression levels. Initially, the expression levels of these hub genes were measured in the dataset GSE12649. Subsequently, two hub genes were validated in the independent datasets GSE5392 and GSE100927. The expression levels of the hub genes were visually represented in boxplots generated using the “ggplot2” package in R. Furthermore, the predictive and discriminatory abilities of the hub genes were evaluated through receiver operating characteristic (ROC) analysis using the “pROC” package, with the area under the curve (AUC) values determined.

Following the identification of hub genes, single-gene gene set enrichment analysis (GSEA) was conducted to uncover their potential functions, utilizing the “clusterProfiler” package. This analysis was carried out using the human reference genome. Subsequently, enrichplot was utilized to visualize the top 5 activating and inhibiting pathways for each gene in the two disease groups.

For immune infiltration analysis, single-sample gene set enrichment analysis (ssGSEA) was carried out using the “GSVA” package to evaluate the relative infiltration levels of 28 immune cell types in each sample of BD. A comparison of immune cell content between the BD and control groups was conducted using the Wilcoxon test. Furthermore, Spearman correlation analysis between infiltrating immune cells and hub genes was performed using the “corrplot” package in R. The resulting correlations between hub genes and immune cells were visualized using lollipop plots.

R software version 4.2.1 was used to perform statistical analyses. Using the Student’s t-test, continuous variables were compared between two groups. A p-value less than 0.05 was statistically significant.

After standardizing the microarray results, the differentially expressed genes (DEGs) were screened by “Limma” package (p < 0.05 and |log FC| > 0.2). The GSE12649 dataset contained 436 DEGs, including 201 upregulated genes and 235 downregulated genes ( Figures 2A, B ). Then, we intersected the resulting data with 2006 ASRGs and a total of 67 differentially expressed arteriosclerosis-related genes (DE-ASRGs) were identified ( Figure 2C ).

In order to analyze the biological functions and pathways, perform GO and KEGG pathway enrichment analysis on the DEGs, ASRG, and DE-ASRG genomes using DAVID ( Figures 2D-F ). Upon comparing the enrichment results of the three genomes, we observed that DE-ASRGs were primarily focused on signal transduction in biological processes. Additionally, negative regulation of apoptotic processes, negative regulation of translation from RNA polymerase II (pol II) promoter, and apoptotic processes were identified as common biological processes across the three genomes. The GO-CC pathways were principally associated with the cytosol, nucleus, and cytoplasm. The GO-MF analysis revealed enrichment in protein binding, identical protein binding and protein homodimerization activity. Additionally, KEGG analysis showed that these DE-ASRGs were enriched in pathways in cancer, PI3K-Akt signaling pathway and human papillomavirus infection (HPV).

To explore the protein interactions among DE-ASRGs, we performed a PPI network containing 67 nodes and 82 edges with medium confidence (score > 0.4) using the STRING database and visualized the network with Cytoscape software, as shown in Figures 3A, B . Subsequently, the MCODE plugin identified a cluster comprising 5 nodes (CSF1R, HIF1A, TIMP3, ITGA7, SPP1) as illustrated in Figure 3C . The cytoHubba plugin was used to score each node gene by 9 randomly selected algorithms and the top 10 hub genes from each algorithm were identified. Five common genes (CSF1R, SPP1, HIF1A, SOX9, and CX3CR1) were selected from the 9 algorithms, as shown in Figure 3D , using the “UpSetR” package. Further screening of common genes through regression analysis resulted in the identification of three candidate hub genes (HIF1A, SOX9 and CX3CR1) as depicted in Figure 3E .

In this study, we conducted WGCNA to identify gene modules associated with BD in the GSE12649 dataset. Initially, the dataset underwent outlier screening, followed by clustering of the remaining samples as depicted in Figure 4A . To establish scale-free networks, a soft-threshold power of β = 6 (yielding a scale-free R² of 0.85) was selected for BD, as shown in Figure 4B . Using a minimum module size of 30, we identified 12 modules, each represented by a unique color ( Figure 4C ). Heat maps illustrating module-trait relationships based on Pearson correlation coefficients were generated to evaluate the connection between each module and BD ( Figure 4D ). Among the 11 modules analyzed, the pink module exhibited the strongest correlation with BD (correlation coefficient = -0.29, P = 0.02), encompassing a total of 186 genes. Subsequently, the 186 genes within the pink module were intersected with the DE-ASRGs, resulting in the identification of 4 candidate hub genes (CTSD, IRF3, NPEPPS and ST6GAL1) ( Figure 4E ). These candidate hub genes were then integrated with those identified through Cytoscape analyses, ultimately yielding a total of 7 candidate hub genes (CTSD, IRF3, NPEPPS, ST6GAL1, HIF1A, SOX9 and CX3CR1).

Lasso regression method, random forest and SVM-RFE method were used to further screen the expression matrix of the 7 candidate hub genes. The LASSO regression model was designed based on BD as well as control samples. By analyzing the Lasso coefficient profiles and selecting the optimal tuning parameter, λ was determined to be 0.0168688 ( Figure 5A ). Subsequently, five candidate genes were identified through the Lasso analysis. Following this, the seven candidate hub genes were fed into the RF classifier, and the top three genes were chosen based on the mean decrease accuracy and mean decrease Gini scale. Three genes were selected as candidate genes from the RF results ( Figure 5B ). Furthermore, SVM-RFE analysis was performed, revealing that the model incorporating three genes exhibited the best Mean Absolute Error (MAE) and Root Mean Square Error (RMSE) ( Figure 5C ). By comparing the overlapping genes obtained from three machine learning methods, two hub genes were selected: CX3CR1 and ST6GAL1 ( Figure 5D ).

We developed SVM models using the BD datasets GSE12649 and GSE5392 to assess the predictive capability of the two hub genes, CX3CR1 and ST6GAL1. The model’s performance on the GSE12649 test set showed an accuracy of 78.90% and a Kappa of 0.58. For the GSE5392 test set, the SVM model achieved an accuracy of 70.83% and a Kappa of 0.42. The SVM results are detailed in Figures 5E, F , as well as in Table 2 .

We investigated the expression of these two genes in the GSE12649 dataset and GSE5392 validation dataset, comparing BD and control samples. As shown in Figure 6A , CX3CR1 and ST6GAL1 showed a significantly lower expression in BD, so as the GSE5392 validation dataset ( Figures 7A, C ). Likewise, we used the AS dataset GSE100927 for validation, and our findings revealed that the expression level of hub genes was higher in the AS group ( Figure 7E ). ROC analysis was used to verify the specificity and sensitivity of the two hub genes for BD and AS diagnosis. The AUC areas of the two hub genes in GSE12649 and GSE5392 were around 0.7, indicating their diagnostic value ( Figures 7B, D ). In GSE100927 dataset, the AUC areas of two hub genes were > 0.8, showing that the hub genes had good diagnostic values.

Single-gene GSEA analysis was performed on these two hub genes to obtain the related pathways of each gene. The results showed its correlation with proteasome, basal cell carcinoma, fatty acid degradation, mineral absorption and so on ( Figures 7G-J ). The interaction network consisted of two hub genes and 23 TFs (Transcription Factors) ( Figure 6A ). ST6GAL1 was regulated by eleven TFs, while CX3CR1 was regulated by fifteen TFs. Moreover, FOXC1, GATA2, and CREB1 were found to interact with both hub genes, indicating a potential close interaction between these TFs and the hub genes.

To further understand the involvement of the hub genes in immune infiltration, we conducted Spearman correlation analysis to investigate the potential relationships between these hub genes and immune cell infiltration. The results showed significantly higher levels of immature B cells, neutrophils, central memory CD4 T cells, and immature dendritic cells, as well as significantly lower levels of activated CD8 T cells in the BD group compared to the control group, as seen in the box plot ( Figure 6B ). Our findings regarding the hub genes, CX3CR1 and ST6GAL1, were consistently exhibited significantly positive correlations with eosinophils, while demonstrating significantly negative correlations with immature B cells ( Figures 6C, D ).