The data analysis process of our study is shown in Figure 1. Three microarray datasets (GSE100927, GSE28829, and GSE208668) were downloaded from the GEO database¹ (Barrett et al., 2013) of the National Center for Biotechnology Information (NCBI) and a single-cell RNA sequencing (scRNA-seq) dataset (GSE253903). The basic information about the datasets is shown in Supplementary Table S1. The GSE100927 dataset contains 69 human arterial samples with AS lesions and 35 without (Steenman et al., 2018). The GSE28829 dataset contains 13 samples of early AS plaques and 16 samples of advanced (Döring et al., 2012). The GSE208668 contains 17 mononuclear cells of peripheral blood samples from patients with ISM and 25 from healthy individuals (Piber et al., 2022). In GSE208668, patients included were aged 60 years and above, with insomnia ≥3 times per week for >3 months, and had not suffered from other sleep disorders or chronic diseases. Removal of batch effects and normalization in microarray datasets was carried out by the “normalizeBetweenArrays” function of the R software “Limma” package (Supplementary Figure S1). The GSE253903 pre-processed by Cellranger (10X Genomics) contains 6 carotid AS plaques from symptomatic patients (Bashore et al., 2024). We searched and downloaded a collection of ISM-related genes in the GeneCards database² (Safran et al., 2010) by the keyword “insomnia” and filtered the top 30% of genes (total 1,952) for subsequent analyses according to the “Relevance score.”

WGCNA is a systems biology technique that reveals gene association patterns in different samples and identifies gene sets that have a significant correlation with phenotypes based on associations between gene sets and between gene sets and phenotypes. After normalization and removing the batch effect of the original datase, the identification was performed using the R software “WGCNA” package. The “goodSamplesGenes” function was used to check the unqualified genes and samples, and then the “pickSoftThreshold” function was used to pick the appropriate soft threshold power (β = 2) to construct the “unsigned” co-expression pattern. Finally, the gene modules were discovered using hierarchical clustering, and a cluster dendrogram was obtained using “plotDendroAndColors.” The module eigengenes (MEs) of different modules were obtained in the first principal component of the modules. Then, the module-trait correlation (Pearson correlation) was assessed based on the association between MEs and the clinical traits. The modules of highest relevance to the trait were selected, Module Membership (MM: the correlation coefficient of gene expression and MEs) and Gene Significance (GS: the correlation coefficient of gene expression and trait) coefficients were calculated (the thresholds: MM > 0.80 and GS > 0.50).

The “Limma” package was used to identify the DEGs in the gene set obtained by WGCNA. The significance threshold was set at adjusted p < 0.05 and | log2 (fold change) | > 0.50. The filter results were visualized in a volcano plot and a heatmap by the “ggplot2” package and “pheatmap” package. The “VennDiagram” package was used to visualize shared genes between gene sets.

The Gene Ontology (GO) term (BP, biological process; CC, cellular component; and MF, molecular function) enrichment analysis (significance p < 0.05 and q < 0.05), the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis (p < 0.05 and q < 0.05), and the Gene Set Enrichment Analysis (GSEA) (p < 0.05 and q < 0.05) were carried out using the “clusterProfiler” package (version 4.10.0). The enrichment results were visualized using the “ggplot2” package and “enrichplot” package.

We constructed the PPI network on the STRING database (version 12.0) (Szklarczyk et al., 2023)³ (Minimum interaction requirement score: medium confidence = 0.400) and imported the results of the network into Cytoscape software (version 3.10.1) (Otasek et al., 2019) for visualization. To filter key sub-network modules from the PPI network, we utilized the Molecular Complex Detection (MCODE) plugin in Cytoscape (parameters: degree cutoff = 2, Node Score Cutoff = 0.2, K-Core = 2, max.Depth = 100). The genes in the key module were identified as the hub genes, which were used to further filter the diagnostic markers.

The least absolute shrinkage and selection operator (LASSO) regression model (a logistic regression method for filtering variables to enhance predictive performance) was constructed based on the “glmnet” package, and the minimum target covariate mean (lambda.min) was determined to filter the candidate biomarkers. The random forest (RF) model was constructed by the “randomForest” package, and the “MeanDecreaseAccuracy” and “MeanDecreaseGini” scores were calculated as the importance scores of the candidate biomarkers. Finally, the shared hub genes screened by the two models were used as candidate biomarkers.

Using the “ggpubr” package, the expression of candidate biomarkers in the GSE208668 and GSE28829 validation sets was compared and visualized in the box diagram. Using the “pROC” package, the diagnostic value of the diagnostic biomarkers was assessed by constructing the ROC, calculating the area under the curve (AUC), and calculating the 95% confidence interval. Nomogram was constructed using the R software “rms” package and its diagnostic efficiency was assessed by calculating the AUC.

Cell-type identification by estimating relative subsets of RNA transcript (CIBERSORT) is a computational method used to translate a normally differentiated gene expression matrix into an infiltrating immune cell proportion. Using the CIBERSORT.R script (Newman et al., 2015), the relative proportions of the infiltration of 22 immune cells in each sample were calculated and depicted in the bar diagram. The comparison of the differential expression of each immune cell between the AS group and the controls is shown by the box diagram. The correlation between diagnostic biomarkers and immune cells (spearman correlation) was calculated using the “corrplot” package.

After downloading the dataset (GSE253903), it was analyzed by the “Seurat” package (version 5.0.1). Quality control of the data was carried: gene counts per cell in the range of 200–2,500 and a percentage of mitochondrial genes less than 5%. Next, the data were normalized using the NormalizeData function, the first 2000 highly variable genes were selected using the “vst” method in the FindVariableFeatures function, the data were scaled using the ScaleData function, and the clustering analysis was performed using the RunPCA function for cluster analysis. We utilized the “harmony” package (version 1.2.0) to eliminate batch effects. Clustering and dimensionality reduction were carried out using FindNeighbors, FindClusters, and RunTSNE functions (dim = 1:20, resolution = 0.5) (Supplementary Figures S2B,C). Subsequently, the cell clusters were visualized using the DimPlot function for visualization. The FindAllMarkers function was utilized to identify the top three marker genes for per cluster. Cell cluster annotation was performed using a strategy of automatic annotation combined with manual correction. Firstly, cell types were annotated using the “SingleR” package (version 2.4.1), and HumanPrimaryCellpronasData downloaded from the “celldex” package (version 1.12.0) was used as the reference datasets. Then, cell types were identified in the CellMarker 2.0 database (Hu et al., 2023)⁴ using the marker genes for each cell cluster. The cell types of each cluster were visualized using DimPlot and DotPlot functions. We utilized the FeaturePlot and VlnPlot functions for visualization to clarify the cellular localization of diagnostic biomarkers.

R software (version 4.3.2) (R Core Team, 2023)⁵ and RStudio software (version 2023.12.1)⁶ were used for data analysis and drawing. The Wilcoxon test was used to compare the differences between the two groups in the test sets. Statistical significance was inferred when p-value < 0.05.

The most strongly related gene modules and gene sets in the AS samples were identified by WGCNA. The most acceptable soft thresholding power (β = 2, R² = 0.930) was chosen based on scale independence and average connectivity (Figure 2A). A total of 7 gene modules were generated in the soft threshold, and the cluster tree diagram of the modules was constructed (Figure 2B). Then, we evaluated the correlation between gene modules and AS (Figure 2C). The results showed the strongest positive correlation between the turquoise module (ME turquoise) and AS (correlation coefficient r = 0.71, p = 3e-17), whereas the blue module (ME blue) had the strongest negative correlation (r = −0.75, p = 1e-19). In addition, we found a strong correlation between MM and GS in the turquoise module (r = 0.71, p < 1e-200) and the blue module (r = 0.81, p < 1e-200) (Figures 2D,E), which indicated a significant correlation between the module genes and AS. Therefore, we filtered the total of 1,502 genes (MM > 0.8 and GS > 0.5) from the two key modules for subsequent analysis.

To filter out genes with a stronger correlation with AS, we analyzed the differential expression of the key module genes. A total of 1,069 DEGs were identified, of which 747 were up-regulated and 322 were down-regulated (Figure 2F). The heatmap demonstrated that the expression of these genes was significantly different in the AS and control groups (Figure 2G). Therefore, the DEGs were identified as a set of AS-related genes.

After taking the intersection of the set of ISM-related genes searched from the GeneCards database with the set of AS-related genes, 61 genes were identified as the shared genes of ISM and AS (Figure 3A). Enrichment analyses of shared genes included GO, KEGG, and GSEA. The GO analysis revealed that the shared genes were majorly enriched for the following functions: (1) biological processes, regulation of innate immune response, regulation of inflammatory response, regulation of neuron projection development, and immune response-activating signaling pathway; (2) cellular component, secretory granule lumen, cytoplasmic vesicle lumen, vesicle lumen, and endocytic vesicle; (3) molecular function, cytokine receptor activity, immune receptor activity, peptide binding, and amide binding (Figures 3B–D). The KEGG analysis revealed that the shared genes were majorly enriched in the following pathways: cytokine-cytokine receptor interaction, efferocytosis, Janus kinase-signal transducer and activator of transcription (JAK–STAT) signaling pathway, chemokine signaling pathway, and lipid and atherosclerosis (Figure 3E). The GSEA revealed that the shared genes were significantly enriched in JAK–STAT signaling pathway and cytokine-cytokine receptor interaction (Figure 3F). These results indicate that inflammation and immune responses may have a vital effect on the development of AS in ISM patients.

To identify the hub genes of ISM-related AS, we uploaded the above shared genes to the STRING database, constructed a PPI network, and visualized it in Cytoscape software. After removing the genes that did not interact with other shared genes, a PPI network containing 41 nodes and 89 edges was constructed (Figure 4A). Then, 4 key modules were identified from the PPI network using the MCODE plugin in Cytoscape (Figures 4B–E). Briefly, module 1 was comprised of 7 genes, including CD4, AIF1, IL10RA, CCR1, HCK, SPI1, and CSF3R. Module 2 comprised 3 genes, including CASP1, APOE, and NCF1. Module 3 comprised 3 genes, including MSH5, XRCC3, and BLM. Module 4 comprised 3 genes, including GNS, GLB1, and HEXB. Obviously, module 1, containing 7 nodes and 21 edges, had the most complex interrelationships. Thus, it was identified as the hub module, and the 7 genes in module 1 were identified as the hub genes. The hub genes have complex interactions with each other, which may have a key role in ISM-accelerated AS.

The hub genes of ISM and AS may contribute to the diagnosis of ISM-related AS patients; therefore, we further selected the diagnostic biomarkers using machine learning. The LASSO regression analysis identified the 6 genes with the lowest binomial deviation among the hub genes with the best fit to the regression model (Figures 5A,B). In addition, we screened the hub genes using the RF algorithm and selected the genes with the top 5 of “MeanDecreaseAccuracy” and “MeanDecreaseGini” scores (Figures 5C,D). After intersecting the results of the above two machine learning algorithms, the 4 candidate diagnostic biomarker genes were selected, including AIF1, IL10RA, CCR1, and SPI1 (Figure 5E).

All 4 candidate diagnostic markers (AIF1, IL10RA, CCR1, and SPI1) were significantly up-regulated in the test set GSE100927 (AS), and the AUCs of the ROCs were over 0.9 (Figures 6A,B). However, in validation sets GSE208668 (ISM) and GSE28829 (AS), it was found that the expression trend of AIF1 was different in the two diseases and the differential expression in GSE28829 was not statistically significant, while IL10RA, CCR1, and SPI1 were up-regulated and the differences were statistically significant (Figures 6C,D). Then, the diagnostic value was assessed by constructing the ROCs and calculating the AUCs. Notably, IL10RA, CCR1, and SPI1 had satisfactory diagnostic value: in GSE208668, IL10RA, AUC = 0.779, 95%CI: 0.639–0.919; CCR1, AUC = 0.786, 95%CI: 0.647–0.925; SPI1, AUC = 0.901, 95%CI: 0.809–0.994; in GSE28829, IL10RA, AUC = 0.918, 95%CI: 0.800–1.000; CCR1, AUC = 0.952, 95%CI: 0.857–1.000; SPI1, AUC = 0.837, 95%CI: 0.670–1.000 (Figures 6E,F). Therefore, we finally identified IL10RA, CCR1, and SPI1 as diagnostic biomarkers. Moreover, to enhance the feasibility of clinical application, the 3 diagnostic markers were utilized to construct the nomogram (Figure 6G). In the nomogram, the expression level of each gene was scored accordingly, and the total score was used to predict the probability of AS. Finally, in the test and validation datasets of AS, the AUC of the nomogram was 0.931 (95% CI: 0.882–0.979) and 0.745 (95% CI: 0.545–0.946), indicating satisfactory diagnostic efficacy (Figure 6H).

Since the functional and pathway enrichment results of the shared genes indicated a potential link between inflammatory-immune with ISM-related AS, we analyzed the characteristics of 22 types of immune cells in AS samples using the CIBERSORT algorithm (Figure 7A). In particular, the proportions of B cells memory, T cells regulatory (Tregs), T cells gamma delta, Macrophages M0, and Mast cells activated were significantly increased in AS samples (all p < 0.05), while the proportions of B cells naive, plasma cells, T cells CD4 memory resting, T cells CD4 memory activated, monocytes, macrophages M1, dendritic cells activated, and mast cells resting were significantly decreased (all p < 0.05) (Figure 7B). Moreover, we analyzed the correlation between gene expression of the three diagnostic biomarkers and the proportions of immune cell infiltration and found that the expression of IL10RA, CCR1, and SPI1 were all significantly positively correlated with the proportions of macrophages M0, T cells gamma delta, and mast cells activated (all p < 0.05) and negatively correlated with the proportions of T cells CD4 memory resting, monocytes, and macrophages M1 (all p < 0.05) (Figure 7C). The diagnostic biomarker genes had the highest correlation coefficients with macrophages M0 proportion (IL10RA, 0.758; CCR1, 0.695; SPI1, 0.774). These results suggest that there may be important interactions between diagnostic biomarkers and various immune cells.

To more precisely describe the immune cell traits and identify the cell types with significant expression of IL10RA, CCR1, and SPI1 in AS plaques, we carried out bioinformatics analysis of the scRNA-seq dataset (GSE253903). After quality control, screening, normalization, and removal of batch effects from the raw data (Supplementary Figure S2A), we performed dimensionality reduction and clustering analysis based on gene expression profiles, and finally obtained a total of 17 cell clusters (Figure 8A). Then, cell annotations identified 11 cell types: common myeloid progenitor (CMP), neutrophil, monocyte, macrophage, myeloid dendritic cell (mDC), plasmacytoid dendritic cell (pDC), B cell, T cell, natural killer (NK) cell, endothelial cell, and smooth muscle cell (Figures 8B,C). Cell expression analysis of the diagnostic marker genes revealed that IL10RA was mainly expressed in monocyte, mDC, pDC, B cell, T cell and NK cell, CCR1 in Monocyte and mDC, and SPI1 in CMP, neutrophil, monocyte, macrophage and mDC (Figures 8D,E). These results indicated that the diagnostic biomarker genes were predominantly expressed in immune cells in AS plaques, with CCR1 and SPI1 expression especially clustered in myeloid immune cells.