The processed data files of microarray were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), a global publicly available repository for genomic data submitted by the scientific community (6). The GSE28829 dataset, tested on GPL570 based on Affymetrix Human GenomeU133 Plus 2.0 Array, is the gene expression data of complete transcriptomics analysis of atherosclerotic carotid artery segments. In total, 13 early lesions (pathological intimal thickening and intimal xanthoma), together with 16 advanced lesions (thin or thick fibrous cap atheroma) are enrolled in this study (7). Gene probes were matched to the official gene symbols after the platform reference matrix files were downloaded. For the situation that more than one probe was matched to a gene symbol, we retained the average gene expression value. The following procedures were analyzed according to the processed matrix file. The flowchart of our analysis was shown in Figure 1.

The normalized gene expression data with the matched gene symbols were utilized to perform Gene Set Enrichment Analysis (GSEA) for interpreting the biological insights of gene expression data (8). We downloaded the GSEA software (http://software.broadinstitute.org/gsea/index.jsp, version 4.0.3) and the reference gene sets C2 (c2.cp.kegg.v7.4.symbols.gmt) from the Molecular Signatures Database (MSigDB, http://www.gsea-msigdb.org/gsea/downloads.jsp) to perform GSEA analysis. The normalized enrichment scores (NES) and nominal p (NP) values were generated automatically by GSEA. |NES|≥1, p-value <0.05, and FDR <0.25 were considered statistically significant. The CIBERSORTx algorithm was applied to infer the relative proportions of infiltrating immune cells of the atherosclerotic samples (9). The matched expression file of GSE28829 was uploaded to the CIBERSORTx online web (http://cibersortx.stanford.edu), with the algorithm run by setting the default signature matrix at 1,000 permutations for quantifying immune cell fractions from bulk tissue gene expression profiles. The CIBERSORTx generates a p-value for the deconvolution for each sample in the results using the Monte Carlo sampling. The significant alteration of immune cells between early and advanced atherosclerosis was identified according to the threshold of the Wilcoxon test at p < 0.05.

The limma package was used to screen DEGs between early and advanced atherosclerosis based on the R platform (R-project.org, version 3.6.1) (10). Fold change (FC) by logarithmic operations with 2 as base numbers were used to make easier calculations and more scientific comparisons for the expression of each gene between early and advanced atherosclerotic samples. Genes with |log₂FC|≥1 were considered as DEGs, and statistical significance was defined by adjusted p-value < 0.05, corrected by the Benjamini-Hochberg method.

Search Tool for the Retrieval of Interacting Genes (STRING, http://string-db.org/, version 11.5) database was used to analyze the protein-protein interaction (PPI) of the DEGs (11). The DEGs were entered into the identifier and Homo sapiens were selected as the organism. To further narrow the candidate gene field, a high confidence level of 0.700 was assessed. Next, the PPI networks were constructed and visualized using the Cytoscape software (version 3.8.2) (12). The plug-in Molecular Complex Detection (MCODE, version 2.0.0), an automated kit to find densely connected regions based on topology, was used to screen the co-expression network of atherosclerosis(13). The MCODE parameters criteria were set by default as follows: degree cut-off = 2, node score cut-off = 0.2, Max depth = 100, and k-score = 2.

The gene ontology (GO) analysis is an international gene functional classification system for determining the correlation between the selected genes and standardized categories through a hypergeometric test (14). The GO enrichment analysis of the cluster genes such as biological process (BP), cellular component (CC), and molecular function (MF) was performed via http://geneontology.org (15). Particularly, Homo sapiens were selected in order to limit the annotation of the species. A p-value of < 0.05 was considered a threshold value.

Pearson correlation analysis was applied to explore the relationship between gene expressions and the relative proportion of immune cells in early and advanced atherosclerosis samples analyzed by the CIBERSORTx (16). Based on the paired t-test, p < 0.05 was considered statistically significant.

Lasso Cox regression analysis was performed via the glmnet package in R software to calculate and select the linear models and preserve valuable variables (17). The expression level of cluster DEGs and the diagnosis of 29 samples were obtained from the probe-matched matrix file, and the samples were randomly assigned to training or testing cohorts in approximately a 2:1 ratio. Variables with calculated coefficients were used to build the classification model. The pROC package in R software was used for the display of Receiver Operating Characteristics (ROC) analysis and the calculation of the area under the curve (AUC) (18). Thus, we investigated the prediction feasibility of hub genes via the AUC values.

For the validation of the identified hub genes, GSE43292 and GSE41571 datasets were downloaded from the GEO database. GSE43292 contains 32 samples of non-atherosclerotic tissue and 32 samples of atheroma plaque, and GSE41571 contains 5 samples of ruptured plaque and 6 samples of stable plaque (19, 20). The expressions of the cluster genes were extracted from the microarray dataset and analyzed by student's t-test. p-value < 0.05 was defined as statistical significance. The ROC curve was constructed to identify the discriminatory power of selected biomarkers correlated to atherosclerosis progression and plaque rupture.

The TRRUST v2 database (www.grnpedia.org/trrust) was applied to analyze the key regulators in transcriptional factors (TFs)-gene interactions to assess the effect of the overlapped TFs on the expression and functional pathways of the cluster genes (21). The cluster genes also served as promising targets for searching for candidate drugs through the Drug-Gene Interaction database (DGIdb, http://www.dgidb.org/, v4.2.0—sha1 afd9f30b) (22).

First, we performed GSEA analysis considering all the gene expressions in GSE28829, not only those above the cutoff in terms of fold change or significance. We noticed that a large number of annotated Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched in advanced atherosclerosis were correlated with immunoinflammatory responses including “NOD_LIKE_RECEPTOR_SIGNALING_PATHWAY” (ES = 0.7229, NP = 0.0021); “B_CELL_RECEPTOR_SIGNALING_PATHWAY” (ES = 0.7219, NP <0.001); “TOLL_LIKE_RECEPTOR_SIGNALING_PATHWAY” (ES = 0.6823, NP = 0.0041); “CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION” (ES = 0.6731, NP <0.001); “CHEMOKINE_SIGNALING_PATHWAY” (ES = 0.6576, NP = 0.0020); and “NATURAL_KILLER_CELL_MEDIATED_CYTOTOXICITY” (ES = 0.6576, NP <0.001). Additionally, other terms including “CELL_ADHESION_MOLECULES_CAMS” (ES = 0.6203, NP = 0.0144), “LEUKOCYTE_TRANSENDOTHELIAL_MIGRATION” (ES = 0.6100, NP = 0.0082), as well as autoimmune-related terms including “SYSTEMIC_LUPUS_ERYTHEMATOSUS”(ES = 0.7604, NP = 0.0061), “AUTOIMMUNE_THYROID_DISEASE” (ES = 0.7269, NP = 0.0021), and “ASTHMA” (ES = 0.7835, NP = 0.0063) were also annotated with significance, indicating the multifaceted mechanisms underlying the progression of atherosclerosis (Figure 2A, see Supplementary Table 1 for details) (3, 4).

According to the results of GSEA analysis, the immune cells may play a causative role in triggering atherosclerotic plaque progression. Therefore, we performed the CIBERSORTx algorithm to comprehensively characterize the infiltration percentages of 19 detectable subpopulations of immune cells in the early and advanced atherosclerotic samples from GSE28829. The relative percentage of each cell in 13 early lesions and 16 advanced lesions is shown in Figure 2B. Moreover, as shown in Figure 2C, the relative proportions of 8 subtypes of the immune cells were significantly different and objectively detectable between early and advanced atherosclerotic samples. Memory B cells (p < 0.01), gamma delta T cells (p < 0.01), M0 macrophages (p < 0.01), and M2 macrophages (p < 0.01) had a higher fraction in advanced atherosclerosis samples, while naïve B cells (p < 0.05), CD8+ T cells (p < 0.01), regulatory T cells (Tregs, p < 0.001), and monocytes (p < 0.01) contained a lower fraction in advanced atherosclerosis samples. The detailed results of the CIBERSORTx analysis and measures of confidence are shown in Supplementary Table 2.

The available numerical expression values of 13 early lesions and 16 advanced lesions from GSE28829 were used to identify the DEGs. Compared with the early atherosclerosis samples, a total of 171 DEGs including 146 up-regulated and 25 down-regulated DEGs were screened in advanced atherosclerosis samples (|log₂FC| ≥ 1 and adjusted p-value < 0.05). The expression data with gene symbols are shown in Supplementary Table 3.

To construct the PPI network of the identified DEGs, the STRING online database and the Cytoscape software were utilized. A total of 171 DEGs were filtered into the PPI network, including 87 nodes and 242 edges (Figure 3A). Based on the high confidence level of 0.70, 84 genes did not fall into the PPI network. According to the node degree> 10 criteria, the 14 hub genes were TYROBP (degree = 25), ITGB2 (degree = 24), CSF1R (degree = 21), CD163 (degree = 17), C1QB (degree = 16), C1QA (degree = 14), CTSS (degree = 13), FCER1G (degree = 13), LY86 (degree = 13),TLR2 (degree = 13), CD86 (degree=12), LAPTM5 (degree = 11), MMP9 (degree = 11), and CXCR4 (degree = 10). The plug-in kit MCODE was conducted to analyze the significant module, and a module with 9 nodes and 24 edges was selected from the PPI network (Figure 3B), showing that the cluster genes were composed of ITGB2, C1QC, LY86, CTSS, C1QA, CSF1R, LAPTM5, VSIG4, and CD163. The descriptions of the cluster genes were listed in Table 1.

To further investigate the biological interpretations of the densely connected genes, GO analysis of the cluster genes was conducted via geneontology.org database, and a total of 118 significant GO terms was annotated (adjusted p-value of < 0.05, see Supplementary Table 4 for details). Consistent with the results of GSEA analysis, the enriched GO terms were also related to immunoinflammatory responses. The leading GO term based in the highest significance were (1) GO:0006954~inflammatory response (adjusted p = 7.09E-6); (2) GO:0045087~innate immune response(adjusted p = 8.75E-6); (3) GO:0006952~defense response (adjusted p = 1.16E-5); (4) GO:0001851~complement component C3b binding (adjusted p = 0.0002); (5) GO:0002526~acute inflammatory response (adjusted p = 0.0005); (6) GO:0043030~regulation of macrophage activation (adjusted p = 0.0006); (7]) GO:0001846~opsonin binding (adjusted p = 0.0008); (8) GO:0002253~activation of immune response (adjusted p = 0.0008); (9) GO:0016485~protein processing(adjusted p = 0.0008); (10) GO:0098883~synapse pruning (adjusted p = 0.0008); (11) GO:0006955~immune response(adjusted p = 0.0008); (12) GO:0042116~macrophage activation (adjusted p = 0.0009). The relations between cluster genes and the leading GO terms are shown in Figure 3C.

The results of previous GSEA analysis and GO enrichment annotation both indicated that immune cells may play an important role in the process of atherosclerosis. Therefore, we carried out a Pearson correlation analysis to illustrate the relationship between cluster genes and immune cells in early and advanced atherosclerotic samples. As shown in Figure 4, it is worth noting that at the threshold of R>0.5 and p < 0.05, almost all the cluster genes represented negative correlations with Tregs and positive correlations with M0 macrophages. Besides, ITGB2, C1QC, LY86, CTSS, C1QA, CSF1R, and CD163 showed negative correlations with CD8+ T cells, similar correlations were also observed in C1QC, CTSS, C1QA, and LAPTM5 with naïve B cells. Moreover, C1QC, LY86, CTSS, and C1QA also exerted positive correlations with M2 macrophages. These results further illustrated that cluster genes may promote atherosclerosis via the regulation of immune cell infiltration and differentiation.

Firstly, the Lasso Cox regression model for the cluster DEGs of atherosclerotic samples from GSE28829 was conducted to investigate an optimum linear combination in predicting advanced atherosclerosis (Figures 5A,B), with coefficients 0.6626 and 0.0228 for ITGB2 and C1QA, respectively. Then, the ROC curve analysis of the Lasso Cox regression classifier was conducted separately to predict advanced atherosclerosis in the training cohort, testing cohort, a combination cohort, and the AUC values all beyond 0.9 (Figure 5C), which suggests it might have a decent potentiality of being biomarkers for screening the progression of the atherosclerosis.

To further validate the accuracy and reliability of the cluster genes selected from GSE28829, two external cohorts from the GSE43292 and GSE41571 datasets containing were included in our study. The detailed information and pathological relationships between the two datasets were shown in Figure 6A, and the expression values of the cluster genes were extracted and analyzed independently. Consistent with the results of GSE28829, these cluster genes (n = 9) were all significantly upregulated in atheroma plaque samples (Figure 6B). We further conducted ROC analysis to investigate the feasibility of C1QA and ITGB2 as prognostic DEGs, and the results showed good predictive accuracy (Figure 6C). For C1QA, the AUC was 83.4961% at the optimal cut-off value of 9.115, and the sensitivity and specificity were 84.375% and 68.75%, respectively; for ITGB2, the AUC was 83.3984% at the optimal cut-off value of 10.03, with the sensitivity being 75% and specificity being 81.25%.

Similarly, in GSE41571, nine (n = 9) of these cluster genes in ruptured plaques were all observed upregulated with seven (n = 7) of them being significant except for VSIG4 and CD163, as compared with stable plaques (Figure 6D). For C1QA, the sensitivity was 100% and specificity was 83.3333% at the optimal cut-off value of 0.8650, with the AUC being 90.6667%; for ITGB2, at the optimal cut-off value of 1.0257, the sensitivity and specificity were 100% and 66.6667%, respectively, and the AUC was 80% (Figure 6E).

Together, the expression levels of cluster genes were validated in two external cohorts, which further demonstrated their importance in the development of atherosclerosis. Additionally, C1QA and ITGB2 could not only differentiate early and advanced atherosclerosis but also highly specific and sensitive DEGs for predicting plaque rupture.

We then attempted to investigate the key regulators of the nine cluster genes, which may be promising targets for atherosclerosis intervention. Three (n = 3) key TFs were identified via the TRRUST v2 database to potentially regulate the expression of CD163, CSF1R, CTSS, and ITGB2 (Table 2). The drug-gene interaction network of the cluster genes was screened aiming to identify druggable targets. As shown in Table 3, a total of 6 intersecting drugs targeting 3 genes, CSF1R, CTSS, and ITGB2, were selected as candidate druggable molecular targets for atherosclerosis.