Raw gene expression profiling of atherosclerosis patients was accessed from GSE28829 (Döring et al., 2012), GSE43292 (Ayari and Bricca, 2013), GSE41571 (Lee et al., 2013), and GSE120521 (Mahmoud et al., 2019) datasets of the Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/gds/). In the GSE28829 dataset, there are 13 early carotid atherosclerotic plaque specimens (pathological intimal thickening and intimal xanthoma) and 16 advanced carotid atherosclerotic plaque specimens (thin or thick fibrous cap atheroma), detected by the Affymetrix Human Genome U133 Plus 2.0 Array. The GSE43292 dataset includes 32 early-stage and 32 advanced-stage carotid atherosclerotic plaque specimens, detected by the Affymetrix Human Gene 1.0 ST Array. The GSE41571 dataset contains five ruptured atherosclerotic plaque specimens and six stable atherosclerotic plaque specimens, detected by the Affymetrix Human Genome U133 Plus 2.0 Array. The GSE120521 dataset comprises four stable and four unstable atherosclerotic plaque specimens. Due to the similarity of sequencing methods, stages of plaque investigated, and study design/comparison between the GSE28829 and GSE43292 datasets, the expression profiling of the aforementioned datasets was merged as the discovery set and batch effects were directly adjusted for batch effects utilizing Combat function of sva package (Leek et al., 2012). Principal component analysis (PCA) was applied for evaluating the performance of the Combat function. The GSE41571 and GSE120521 datasets were utilized as the external verification sets. The probe ID for each gene was transformed into a gene symbol. If a gene symbol corresponded to several probe IDs, the average expression value of the probe IDs was calculated as the representative expression value of the gene.

The list of 1,242 immune-associated genes was curated from the Immunology Database and Analysis Portal (ImmPort; https://www.immport.org/home) (Bhattacharya et al., 2018). Through limma package (Ritchie et al., 2015), differentially expressed immune-associated genes were screened with 45 early-stage and 48 advanced-stage carotid atherosclerotic plaques in line with the criteria of |fold-change|>1.5 and false discovery rate (FDR) < 0.05. These genes were regarded as atherosclerotic plaque progression- and immune-associated genes.

ClusterProfiler package (Yu et al., 2012) was utilized for functionally analyzing the biological functions, which comprises Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). The p-value was adjusted using the Benjamini–Hochberg approach or FDR for multiple testing corrections. The threshold was set at FDR<0.05. GO categories comprised biological processes (BP), molecular functions (MF), and cellular components (CC).

The atherosclerotic plaque progression- and immune-associated genes were uploaded onto the online “Search Tool for the Retrieval of Interacting Genes” (STRING; http://string-db.org) and their interaction pairs were required. Through the plug-in of Cytoscape Molecular Complex Detection (MCODE) (Bader and Hogue, 2003), hub modules of the PPI network were established following the threshold of degree cutoff = 2, K-Core = 2, and node score cutoff = 0.2.

Three machine learning algorithms, LASSO, random forest, and SVM-RFE (Sanz et al., 2018), were applied for screening characteristic genes. LASSO, a dimension reduction approach, shows superiority in evaluating high-dimensional data in comparison to regression analysis. LASSO analysis was implemented with a turning/penalty parameter utilizing a 10-fold cross-verification via glmnet package (Engebretsen and Bohlin, 2019). Recursive feature elimination (RFE) from the random forest algorithm, an approach of supervised machine learning, was applied for ranking the atherosclerotic plaque progression- and immune-associated genes. The predictive performance was estimated via ten-fold cross-validation, and the genes with relative importance>0.25 were determined as the characteristic genes. SVM-RFE is superior to linear discriminant analysis and to the mean squared error method to select relevant characteristics and remove redundant characteristics. SVM-RFE was applied for feature selection via ten-fold cross-validation. Receiver operating characteristic (ROC) curves and the area under the curve (AUC) were used for estimating the diagnostic efficacy.

Single-sample gene set enrichment analysis (ssGSEA) was implemented to analyze the infiltration levels of immune cells on the basis of the expression profiling of 29 immunity-relevant signatures.

GSEA was implemented for functionally elucidating the biological significance of characteristic genes (Subramanian et al., 2005). The gene set of “c2.cp.kegg.v11.0.symbols” from the Molecular Signature Database (MSigDB, http://software.broadinstitute.org/gsea/msigdb) (Liberzon et al., 2015) was utilized as the reference set. For achieving a normalized enrichment score for each analysis, gene set permutations with 1,000 times were conducted. An FDR<0.05 was regarded as significant enrichment.

Characteristic genes were incorporated to establish a nomogram using the rms package. The calibration curve was utilized for evaluating the accuracy of the nomogram. Through the decision curve analysis, the clinical usefulness of the nomogram was evaluated.

The Connectivity Map (CMap, https://clue.io/), a web-based database, applies cellular responses to perturbations for finding interactions between diseases, genes, and small-molecule compounds (Subramanian et al., 2017). The atherosclerotic plaque progression- and immune-associated genes were interrogated to compare the similarity to all perturbed signatures in this database. The candidate small-molecule compounds were determined with an |enrichment score|>90. Moreover, compounds with positive or negative enrichment scores were selected for predicting the mode of action (MoA).

The consensus clustering approach was applied to quantitatively estimate the number of unsupervised classes across carotid atherosclerotic plaque specimens via the ConsensusClusterPlus package (50 iterations and resampling rate of 80%) on the basis of expression profiling of atherosclerotic plaque progression- and immune-associated genes (Wilkerson and Hayes, 2010). The consensus matrix plot, consensus cumulative distribution function (CDF) plot, relative alterations in area under the CDF curve, and tracking plot were implemented for finding the optimal number of clusters. Principal component analysis (PCA) was utilized for defining the expression difference in atherosclerotic plaque progression- and immune-associated genes between two subtypes. The PCA diagram was depicted utilizing the ggplot2 package (Ito and Murphy, 2013).

GSVA is a non-parametric and unsupervised gene set enrichment approach, which evaluates the association between biological pathways and gene signatures on the basis of expression profiling (Hänzelmann et al., 2013). Fifty hallmark gene sets were curated from the MSigDB as the reference set. The GSVA package and its ssGSEA function were applied for obtaining the GSVA score of each gene set. The GSVA score denoted the degree of absolute enrichment of each gene set. Limma package was utilized for comparing the difference in the GSVA score of each gene set between subtypes.

All statistical tests were implemented utilizing R software 3.6.1. Wilcoxon or Student’s t-test was utilized for analyzing the difference between the two groups. The correlation between the variables was determined using Pearson’s or Spearman’s correlation test. All statistical p-values were two-sided, and p < 0.05 was regarded as statistical significance.

To investigate the roles of immune-associated genes in the progression of atherosclerotic plaques, we combined the expression profiles of 45 early-stage and 48 advanced-stage atherosclerotic plaque specimens from the GSE28829 and GSE43292 cohorts (Figure 1A). Batch effects were adjusted for subsequent analysis (Figure 1B). Among 1,242 immune-associated genes, 114 presented downregulation and 21 presented upregulation in an advanced-stage compared to early-stage atherosclerotic plaques (Figures 1C,D). The detailed information is listed in Supplementary Table 1. These atherosclerotic plaque progression- and immune-associated genes were linked to immune responses such as cytokine–cytokine receptor interaction and chemokine signaling pathway (Figures 1E–H).

Through MCODE analysis, one hub module from the PPI network was established, which comprised key atherosclerotic plaque progression- and immune-associated genes (Figure 2A). Further analysis displayed that they mainly participated in cytokine–cytokine receptor interaction, Toll-like receptor signaling pathway, chemokine signaling pathway, TNF signaling pathway, natural killer cell-mediated cytotoxicity, NF-kappa B signaling pathway, IL-17 signaling pathway, and cell adhesion molecules (Figures 2B–E), indicating the crucial roles of key atherosclerotic plaque progression- and immune-associated genes.

Three algorithms were applied for selecting characteristic genes among key atherosclerotic plaque progression- and immune-associated genes. For the LASSO algorithm, the optimal lambda was 0.014 following ten-fold cross-validation. Thus, we chose the minimum criteria for building the LASSO classifier due to higher accuracy by comparisons, and 12 characteristic genes were identified, containing ITGAL, IL7R, IL18, CCL19, BTK, TLR8, CD14, CCL5, IL1B, IL6, GZMB, and TNFSF13B (Figures 3A,B). For the random forest algorithm, 30 characteristic genes with relative importance >0.25 were determined, including MMP9, ICAM1, PTPRC, LCP2, C3AR1, CCL2, IL10RA, IL6, FCGR3A, CD28, TNFSF13B, TLR2, CCL5, CD4, CD86, TLR1, CSF2RB, TYROBP, CCL4, ITGB2, FCER1G, CSF1R, CYBB, CCL19, HCK, CCR1, ITGAL, CD14, GZMB, and BTK (Figures 3C,D). For the SVM-RFE algorithm, when the feature number was 26, the classifier had the minimum error, containing CD14, ITGAL, TNFSF13B, IL18, CCL5, PTPRC, CCRL2, IL7R, MMP9, BTK, IL10RA, CD28, GZMB, ICAM1, HCK, CSF2RB, CD74, TLR2, CCR1, C3AR1, CCL19, IL2RG, TYROBP, CSF1R, CXCR4, and IL1B (Figure 3E). Following intersection, 7 characteristic genes shared by LASSO, random forest, and SVM-RFE algorithms were finally identified (TNFSF13B, CCL5, CCL19, ITGAL, CD14, GZMB, and BTK; Figure 3F).

Seven characteristic genes (TNFSF13B, CCL5, CCL19, ITGAL, CD14, GZMB, and BTK) presented higher expression in advanced-stage than in early-stage atherosclerotic plaques (Figure 3G), indicating their potential roles during the progression of atherosclerotic plaques. When all of them were fitted into one variable, the AUC of the ROC curve was 0.918, demonstrating the favorable diagnostic efficiency in predicting atherosclerotic plaque progression (Figure 3H). We also estimated the diagnostic performance of each characteristic gene in predicting atherosclerotic plaque progression in the combined GSE28829 and GSE43292 cohorts. The AUC values of ROC curves were 0.888 of BTK (Figure 4A), 0.753 of CCL5 (Figure 4B), 0.834 of CCL19 (Figure 4C), 0.854 of CD14 (Figure 4D), 0.862 of GZMB (Figure 4E), 0.872 of ITGAL (Figure 4F), and 0.816 of TNFSF13B (Figure 4G), demonstrating that these characteristic genes enabled to estimate the progression of atherosclerotic plaques.

The expression of characteristic genes was verified in external datasets. In the GSE41571 dataset, CCL19 presented higher expression in stable than in ruptured plaques, while ITGAL, CD14, and GZMB had higher expression in ruptured than in stable plaques (Figure 4H). In the GSE120521 dataset, higher TNFSF13B, CD14, and BTK expression was confirmed in unstable than in stable plaques (Figure 4I). The AUC values of ROC curves were 0.717 of BTK (Supplementary Figure 1A), 0.567 of CCL5 (Supplementary Figure 1B), 0.917 of CCL19 (Supplementary Figure 1C), 0.900 of CD14 (Supplementary Figure 1D), 0.933 of GZMB (Supplementary Figure 1E), 0.900 of ITGAL (Supplementary Figure 1F), and 0.700 of TNFSF13B (Supplementary Figure 1G) in the GSE41571 dataset, indicating their potential in distinguishing ruptured from stable plaques. Moreover, the AUC values of the ROC curves were 1.000 of BTK (Supplementary Figure 2A), 0.812 of CCL5 (Supplementary Figure 2B), 0.938 of CCL19 (Supplementary Figure 2C), 1.000 of CD14 (Supplementary Figure 2D), 0.875 of GZMB (Supplementary Figure 2E), 0.812 of ITGAL (Supplementary Figure 2F), and 0.969 of TNFSF13B (Supplementary Figure 2G) in the GSE120521 dataset, demonstrating that they are capable of differentiating unstable from stable plaques. Hence, the characteristic genes possessed excellent diagnostic performance in predicting the progression of atherosclerotic plaques.

Immunological features were evaluated in accordance with immune cell infiltration and immune checkpoint expression. Compared with early-stage atherosclerotic plaques, most innate and adaptive immune cells presented higher infiltration levels in advanced-stage atherosclerotic plaques (Figure 5A). Moreover, there were remarkable interactions between immune cell populations across atherosclerotic plaques (Figure 5B). As illustrated in Figure 5C, the higher expression of most immune checkpoints was investigated in advanced-stage than in early-stage atherosclerotic plaques. The aforementioned data indicated a higher immune response in advanced-stage atherosclerotic plaques. Furthermore, analyses displayed positive interactions between characteristic genes and immune cell infiltrations (Figure 5D). Additionally, characteristic genes were positively linked to immune checkpoints across atherosclerotic plaques (Figure 5E). Hence, the characteristic genes might modulate immunological features during atherosclerotic plaque progression.

Through GSEA, we evaluated signaling pathways involved in the characteristic genes. Our results demonstrated that BTK (Figure 6A), CCL5 (Figure 6B), CCL19 (Figure 6C), CD14 (Figure 6D), GZMB (Figure 6E), ITGAL (Figure 6F), and TNFSF13B (Figure 6G) were all positively linked to the immune responses (cytokine–cytokine receptor interaction, Toll-like receptor signaling pathway, B-cell or T-cell receptor signaling pathway, etc.).

As illustrated in Figure 7A, there were remarkable interactions between the characteristic genes. By incorporating characteristic genes, a nomogram was constructed as a diagnostic tool for atherosclerotic plaque progression (Figure 7B). In the nomogram, each characteristic gene corresponded to a score, and the total score was calculated by adding the scores for all characteristic genes. The total points corresponded to different risks of atherosclerosis. The calibration curve demonstrated that the nomogram enabled an accurate estimation of the progression of atherosclerotic plaques (Figure 7C). As depicted in the decision curve analysis, the patients diagnosed with atherosclerosis could clinically benefit from the nomogram (Figure 7D).

On the basis of atherosclerotic plaque progression- and immune-associated genes, potential small molecular compounds against atherosclerosis were predicted through CMap analysis, as depicted in Figure 7E. Among them, alvespimycin, pifithrin-mu, and radicicol shared HSP inhibitors, while IKK-2-inhibitor-V and radicicol shared NF-kappa B pathway inhibitors.

Through the consensus clustering approach, atherosclerotic plaques were clustered in accordance with expression profiling of 135 atherosclerotic plaque progression- and immune-associated genes. The optimal number of subtypes was 2, which was determined using a consensus matrix plot, a CDF plot, relative alterations in the area under the CDF curve, and a tracking plot (Figures 8A–D). The two immune subtypes were named C1 and C2. PCA demonstrated the remarkable difference between the subtypes (Figure 8E). It was found that there was remarkable heterogeneity in the expression of atherosclerotic plaque progression- and immune-associated genes between subtypes (Figure 8F).

In Figure 9A, we noticed that all characteristic genes presented a higher expression in C1 than C2 subtype. In comparison to the C2 subtype, most immune checkpoints were remarkably upregulated in the C1 subtype (Figure 9B). As illustrated in Figure 9C, the C1 subtype had higher immune activation (allograft rejection, complement, interferon-gamma response, IL6-JAK-STAT3 signaling, inflammatory response, TNFα signaling via NF-kappa B, etc.) than the C2 subtype. Further analysis demonstrated that the C1 subtype presented higher infiltration levels of most immune cell populations than the C2 subtype (Figure 9D). Collectively, we identified C1 as an immune subtype and C2 as a non-immune subtype.