A total of four AS-related gene expression matrices (GSE100927, GSE43292, GSE28829, and GSE163154) were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). In the GSE100927 dataset (Steenman et al., 2018), atherosclerotic plaques and control arteries without such lesions (obtained from deceased organ donors) were collected from the carotid, femoral, and infra-popliteal arteries. Total RNA was extracted from the samples using standard protocols and subsequently hybridized to microarrays for further analysis. The “limma (Ritchie et al., 2015)” R package was employed to normalize the data and identify differentially expressed genes (DEGs). The significance criteria were defined as −0.58 < log2 fold change (FC) < 0.58, accompanied by a p-value of less than 0.05. In R version 4.2.1, the “ggplot2” and “ComplexHeatmap (Gu et al., 2016)” packages were utilized to generate volcano plots and heatmaps for the differentially expressed genes (DEGs). A list of 332 lactylation-related genes (LRGs) was obtained from previously published reports (Cheng et al., 2023). The collected LRG data were then integrated with the results of differential expression analysis of genes (DEGs) between the disease and control groups in the GSE100927 dataset, and these genes were classified as AS-related LRGs (AS-LRGs). The results were visualized utilizing the “ggplot2” and “VennDiagram” packages.

The GSE43292 dataset was utilized as the validation dataset. In GSE43292 dataset (Ayari and Bricca, 2013), a total of 34 patients who underwent carotid endarterectomy at the University Hospital of Lyon participated in this study. Carotid endarterectomy samples were obtained in the operating room and immediately divided into two fragments: the atheroma plaque and macroscopically intact tissue. We evaluated the expression levels of AS-LRGs in both affected and control groups, using the Wilcoxon rank-sum test for statistical analysis. The ROC analysis on the dataset was performed using the “pROC (Robin et al., 2011)” package. Statistical analysis of paired samples was conducted using a paired t-test. All visualizations were generated using the “ggplot2” package.

Functional and pathway enrichment analyses were conducted using the Gene Ontology (2015) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2022) databases, employing the R package “clusterProfiler (Yu et al., 2012)” based on significant AS-LRGs (P < 0.05).

Gene sets redefined from the MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb/collections.jsp) were utilized to analyze the distribution patterns of genes within a gene expression profile. This profile is ranked based on the correlation between gene expression and phenotypic characteristics, thereby evaluating the contributions and associations of these genes with the phenotype (Subramanian et al., 2005). We retrieved the gene sets “c2.cp.kegg.v2022.1.Hs.symbols” and “c5.all.v2022.1.Hs.symbols” from the MSigDB database and utilized the “clusterProfiler (Yu et al., 2012)” package in R to perform the analysis.

The TarBase v9.0 (Skoufos et al., 2024) and ENCODE (2012) databases were employed to identify the miRNAs and transcription factors associated with AS-LRGs, which were subsequently visualized using the NetworkAnalyst (2012) platform. We utilized the Drug Gene Interaction Database (DGIdb) (Cannon et al., 2024) (https://www.dgidb.org) to predict the relationships between genes and drugs, which were visualized using Cytoscape (Shannon et al., 2003).

The immune infiltration within the dataset was assessed using the single-sample GSEA (ssGSEA) algorithm from the R package GSVA (Hänzelmann et al., 2013), incorporating 24 immune cell markers detailed in an Immunity article (Bindea et al., 2013). We conducted the Wilcoxon rank-sum test to examine the differences between the normal group and the AS group. Spearman’s correlation analysis was conducted to assess the relationship between AS-LRGs and the level of immune infiltration in AS tissues. All visualizations were generated utilizing “ggplot2” package.

In the GSE28829 database, Seventeen samples of atherosclerotic tissue were collected from carotid artery segments, including nine samples representing early lesions and eight samples representing advanced lesions, sourced from the Maastricht Pathology Tissue Collection for inclusion in this study (Döring et al., 2012). The GSE28829 database was utilized to analyze the differential expression of AS-LRGs between early and advanced lesions. In GSE163154 database, they conducted an analysis of the genetic profiles of atherosclerotic lesion segments categorized as low-risk or high-risk in patients undergoing carotid endarterectomy, based on the presence or absence of intraplaque hemorrhage (Jin et al., 2021). The GSE163154 database was employed to analyze the differential expression of AS-LRGs in the absence and presence of intraplaque hemorrhage. The statistical and presentation methodologies utilized in this study are consistent with those previously described.

Figure 1 provides a comprehensive overview of the workflow. The differentially expressed genes (DEGs) in the GSE100927 dataset were visualized using volcano plots (Figure 2A). As observed, there were 1,201 upregulated DEGs and 765 downregulated DEGs. A classification heatmap was created to clearly present these results (Figure 2B). Additionally, by integrating the LRGs with the upregulated and downregulated DEGs from the GSE100927 dataset, we identified 9 upregulated genes (FABP5, LSP1, IKZF1, MNDA, CRABP2, ARID3A, RCC2, WAS, and HMGA1) and 8 downregulated genes (CSRP2, PPP1CB, CSRP1, HIST1H1C, HEXIM1, CALD1, PDLIM1, and RANBP2) that are common to both DEGs and LRGs and are associated with atherosclerosis (Figures 2C, D). Subsequently, we assessed the expression levels of these genes in the GSE43292 dataset. Among the identified genes, FABP5 and HIST1H1C were not detected as independent samples in the GSE43292 dataset. Subsequently, the expression of the 8 upregulated genes and 7 downregulated genes in the GSE43292 dataset were illustrated using box plots (Figures 3A, B). In the validation dataset (GSE43292), the 5 upregulated genes identified as having diagnostic value were: LSP1, IKZF1, MNDA, RCC2, and WAS (AUC >0.7) (Figures 3C, D); the 7 downregulated genes identified as having diagnostic value were: CSRP2, PPP1CB, CSRP1, HEXIM1, CALD1, PDLIM1 and RANBP2 (AUC >0.7) (Figures 3E, F). Subsequently, we conducted paired-sample t-tests on the 5 upregulated genes and 7 downregulated genes within the validation dataset, resulting in significant differences for all genes (Figures 4A–L).

The gene ontology (GO) analysis revealed significant enrichment in a multitude of biological processes, cellular components, and molecular functions. The KEGG analysis of the AS-LRGs revealed significant enrichment across a diverse array of pathways. The results of the GO and KEGG analyses for AS-LRGs are presented in Figure 5A. To evaluate the impact of gene expression levels on AS, Gene Set Enrichment Analysis (GSEA) was conducted using the GSE100927 dataset to elucidate the relationships between gene expression and related biological processes, cellular components, and molecular functions. The findings indicated that in the GSE100927 dataset, the genes were primarily associated with contractile fiber, I band, adaptive immune response, leukocyte-mediated immunity, immune response-regulating signaling pathways, and other biological functions (Figure 5B). The biological pathways predominantly regulated by genes in the GSE100927 dataset include those associated with lysosome, allograft rejection, antigen processing and presentation, leishmania infection, autoimmune thyroid disease, natural killer cell mediated cytotoxicity, cell adhesion molecules, chemokine signaling pathway, cytokine-cytokine receptor interaction, graft versus host disease, asthma, and type I diabetes mellitus (Figures 5C–H).

We constructed an interaction network between AS-LRGs and miRNAs, encompassing 12 genes and 561 miRNAs (Figure 6A). The top five predictions of the AS-LRG-miRNA interaction network included: RANBP2 (targeted by 294 miRNAs), PPP1CB (targeted by 276 miRNAs), CALD1 (targeted by 170 miRNAs), RCC2 (targeted by 164 miRNAs) and HEXIM1 (targeted by 146 miRNAs). The interaction network between AS-LRGs and transcription factors (TFs) consisted of 10 genes and 212 transcription factors (Figure 6B). The top five predictions of the AS-LRG-TF interaction network encompassed: HEXIM1 (targeted by 176 TFs), PPP1CP (targeted by 28 TFs), CSRP1 (targeted by 21 TFs), CALD1 (targeted by 19 TFs) and WAS (targeted by 16 TFs). The AS-LRG–drug interaction network consisted of 11 drugs and 4 genes, as illustrated in Figure 6C.

To achieve a more comprehensive understanding of the role of immune cells in AS and their association with critical genes, we employed the single-sample GSEA methodology to assess immune infiltration within the GSE100927 dataset. Among the 24 immune cell subpopulations, we observed significant differences between AS-affected tissues and normal tissues. The violin plot illustrated that AS patients exhibited varying levels of the majority of immune cell types (Figure 7A). Five upregulated genes and seven downregulated genes were identified as being associated with each immune cell type, as illustrated in Figures 7B, C. The chart indicates that the majority of upregulated genes exhibited a positive correlation with most immune cell types, with the exception of NK cells. Conversely, the chart shows that the majority of downregulated genes exhibited a negative correlation with most immune cell types, with the exception of NK cells.

To evaluate the differential expression of AS-LRGs in early and advanced lesions, we utilized the GSE28829 dataset, which consists of 9 early lesion samples and 8 advanced lesion samples. The expression levels of the five upregulated genes and seven downregulated genes in the GSE28829 dataset were illustrated using box plots (Figures 8A, B). In the GSE28829 dataset, the five upregulated genes identified as having diagnostic value were LSP1, IKZF1, MNDA, RCC2, and WAS (AUC >0.7) (Figure 8C). The five downregulated genes identified as having diagnostic value were CSRP2, PPP1CB, CSRP1, HEXIM1, and CALD1 (AUC >0.7) (Figure 8D).

The GSE163154 dataset, which includes samples indicating the presence or absence of intraplaque hemorrhage, was utilized to assess the relationship between AS-LRGs and intraplaque hemorrhage. The expression levels of the five upregulated genes and seven downregulated genes in the GSE163154 dataset were illustrated using box plots (Figures 9A, B). In the GSE163154 dataset, the five upregulated genes identified as having diagnostic value were LSP1, IKZF1, MNDA, RCC2, and WAS (AUC >0.7) (Figures 9C, D). The six downregulated genes identified as having diagnostic value were CSRP2, PPP1CB, CSRP1, HEXIM1, CALD1 and PDLIM1 (AUC >0.7) (Figures 9E, F).