The entire study process is depicted in Figure 1 . Five gene expression profiles (GSE100927, GSE43292, GSE28829, GSE57691, and GSE159677) related to AS were retrieved from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database using the keywords “atherosclerosis” with species limited to “Homo sapiens.” The GSE100927 dataset (GPL17077) includes 35 healthy control arteries and 69 atherosclerotic diseased arteries, totaling 104 samples that served as the training set (7). The datasets GSE43292 (GPL6244), GSE28829 (GPL570), and GSE57691 (GPL10558) were used as validation sets, with the GSE43292 dataset including 32 healthy control arteries and 32 atherosclerotic diseased arteries, totaling 64 samples (8); the GSE28829 dataset including 13 early atherosclerotic plaques and 16 advanced atherosclerotic plaques, totaling 29 samples (9); and the GSE57691 dataset including 10 control samples and 9 abdominal aortic samples from patients with abdominal aortic occlusive disease, totaling 19 samples (10). Lastly, the GSE159677 dataset (GPL18573) for single-cell data analysis includes 3 carotid atherosclerotic cores and 3 control samples (11).

The GSE100927 dataset was used to identify the DEGs associated with AS by comparing the disease and control groups using the “limma” R package, with the criteria set at adj.P<0.05 and |logFC|>1 for DEGs selection (12). Lipid metabolism-related genes (LMRGs) were obtained from the Molecular Signature Database (MsigDB, https://www.gsea-msigdb.org/) and supplemented with genes having a “relevance score” greater than the mean (4.1689) from the GeneCards (https://www.genecards.org/, accessed on March 15, 2025) database under the keyword “lipid metabolism”. Ferroptosis-related genes (FeRGs), including driver, suppressor, and marker genes, were obtained from the FerrDb V2 (http://www.zhounan.org/ferrdb/current/, accessed on March 15, 2025) database (13). The intersection of DEGs, LMRGs, and FeRGs yielded ASLMRFeGs.

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes (KEGG) were performed using the DAVID (https://davidbioinformatics.nih.gov/home.jsp, accessed on March 15, 2025) database (P<0.05). Gene Set Enrichment Analysis (GSEA) was conducted via the “GSEA” and “clusterProfiler” R package (14). KEGG pathway gene sets (c2.cp.kegg.v7.4.symbols.gmt) were obtained from the MSigDB database. Genes were ranked in descending order based on logFC values of differentially expressed genes. GSEA was performed using the clusterProfiler R package with default weighted enrichment statistic calculation. Multiple hypothesis testing correction was applied using the Benjamini-Hochberg method, retaining significant pathways with FDR < 0.05. The most significantly enriched pathways (NES > 0 for upregulated, NES < 0 for downregulated) were visualized using enrichplot.

We performed immune cell quantification using the CIBERSORT deconvolution algorithm with the LM22 signature matrix (containing gene expression signatures of 22 human immune cell subtypes) (15). Data preprocessing involved automatic logarithmic transformation, quantile normalization, and Z-score normalization. The ν-support vector regression (ν-SVR) algorithm was implemented with parameter optimization (ν=0.25-0.75), model selection by root-mean-square error minimization, and non-negative constraint enforcement. To evaluate the reliability of our results, we performed a 1000 permutation test (perm=1000), retaining significant results with a P < 0.05. Differences between the two groups (C1 compared to C2 cluster, control compared to AS group) were compared using the Wilcox test, and the results were visualized using the “vioplot” R package (16). Subsequently, Spearman correlation analysis was employed to reveal the relationship between ASLMRFeGs and immune cells.

Unsupervised clustering analysis of AS samples based on ASLMRFeGs expression profiles was performed using the “ConsensusClusterPlus” R package (17). The AS samples were grouped by applying the k-means algorithm with 1000 iterations, k=9, reps=50, pItem=0.8, pFeature=1, clusterAlg=km, distance=euclidean. The appropriate number of clusters was determined based on the matrix heat map, cumulative distribution function (CDF) curve, consensus matrix, and consistent cluster score (>0.9).

We performed pathway activity quantification using GSVA on normalized expression profiles with the “c2.cp.kegg.symbols.gmt” gene set (18). The expression matrix was de-duplicated using the avereps function and filtered to include only AS group samples. Single-sample pathway enrichment scores were calculated using the GSVA algorithm, normalized by min-max scaling. Differential pathway analysis between predefined clusters (C1 compared to C2 cluster) was performed using two-sample t-tests (adj.P<0.05). The final visualization highlighted the up- and down-regulated pathways with the most significant differences (|t-value| > 2).

We performed co-expression network analysis using the WGCNA R package (19). We retained genes with expression standard deviations ranking in the top 25% of highly variable genes and excluded low-expression genes and low-quality samples. The Euclidean distance matrix of gene expression profiles was calculated, and a hierarchical clustering dendrogram was constructed using the average linkage method. A cut height of 80 was set to remove samples that were too distant from the central cluster. The pickSoftThreshold function was used to evaluate the scale-free topology fit, selecting the minimum β value (soft thresholding power) that ensured R² > 0.9 for scale-free network construction. A weighted adjacency matrix was constructed based on the Pearson correlation coefficient raised to the power of the chosen β, and the topological overlap matrix (TOM) was further calculated to reduce noise interference. Dynamic tree cutting clustering was performed with parameters set as follows: minModuleSize=100, deepSplit=2, mergeCutHeight=0.25, and modules were divided using the cutreeDynamic function, with highly correlated modules (cor > 0.75) merged using the mergeCloseModules function. The first principal component of each module was calculated using moduleEigengenes. Pearson correlation with phenotypic data was computed, and significance was adjusted using Student’s t-test. Genes with expression correlated with the target trait (GS) and module eigengenes (MM) were identified, with genes defined as hub genes if GS > 0.5 and MM > 0.8.

To further narrow down the Hub genes, the module Hub genes were imported into the STRING database (https://www.string-db.org/, accessed on March 15, 2025), and “Homo sapiens” was selected as the species with high confidence (0.7) (20). The PPI network of modular Hub genes was constructed, and topological analysis was performed using the Cytoscape plugin “CytoHubba” to filter out the candidate Hub genes with degree centrality, betweenness centrality, and closeness centrality greater than the mean and construct the PPI network of candidate Hub genes (21).

The Least Absolute Shrinkage and Selection Operator (LASSO) regression algorithm is a linear model that achieves feature selection and sparsity by introducing an L1 regularization term, effectively screening out important variables and reducing model complexity (22). Support Vector Machine-Recursive Feature Elimination (SVM-RFE) optimizes the feature set by recursively training SVM models and eliminating the features with the smallest weights (23). Random Forest (RF) is an ensemble learning method that assesses the importance and contribution of features by constructing multiple decision trees.

We utilized the “glmnet” R package to perform LASSO regression with a binomial distribution (family=binomial) and L1 regularization (alpha=1). The optimal lambda value was determined through 10-fold cross-validation using deviance as the evaluation metric. Subsequently, we extracted the non-zero coefficient genes to establish the final feature gene set, excluding the intercept term.

For SVM-RFE analysis, we employed the “e1071,” “kernlab,” and “caret” R packages. A radial basis function kernel was used, and classification performance was evaluated using 10-fold cross-validation with root mean square error (RMSE) as the metric. Feature subset evaluation included specific gene numbers (2, 4, 6, and 8 genes) and continuous sequences (10–40 genes with a step size of 3). The feature gene combination that minimized cross-validation RMSE was selected as the optimal feature set.

RF analysis was conducted using the “randomForest” R package with parameters set as follows: 500 initial decision trees (ntree=500), default mtry parameter (square root of features), and classification mode. The optimal number of trees was determined by minimizing the out-of-bag (OOB) error, followed by model reconstruction. Feature importance was assessed by calculating the mean decrease in Gini index, and genes with importance scores >2 were identified as significant features. The area under the receiver operating characteristic (ROC) curve was visualized using the “pROC” R package (24).

A nomogram model was established using the “rms” R package to predict the risk of developing AS, and its predictive power was estimated by using calibration curves and DCA.

The GSE28829, GES43292 and GSE57691 datasets were used for external validation, and the Wilcoxon test was used to compare the differences in expression between the two groups, with P<0.05 considered significant.

Single-cell sequencing of scRNA-seq data from the GSE159677 dataset was performed using the “Seurat” R package to detect marker genes in each cell cluster (25). Using the “Seurat” R package, initial quality control was performed by excluding cells with fewer than 300 or more than 7500 detected genes, as well as cells with granulocyte gene proportions greater than 15% or hemoglobin gene proportions greater than 1%. Data were normalized using LogNormalize (scale factor=10,000) and the top 2,500 highly variable genes were selected via variance-stabilizing transformation. Dimensionality reduction employed principal component analysis (20 PCs selected by elbow method) followed by cell clustering (Shared Nearest Neighbor algorithm, resolution=0.3) and UMAP visualization (n.neighbors=30, min.dist=0.3). Cell populations were annotated using established marker genes ( Supplementary Table S1 ), and cellular proportions across experimental groups were quantified and visualized using stacked bar plots.

Imported the ASLMRFeGs and Hub genes into the DGIbd database (https://dgidb.org/) to obtain drug candidates and screen for FDA-approved drugs (26).

The 3D structure of drugs was obtained from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/), and the crystal structure of protein was obtained from the PDB database (https://www.rcsb.org/), and the protein receptor and small molecule ligand were dehydrated and hydrogenated using PyMOL software, and then the molecular docking analysis was performed using AutoDockTools and PyMOL software for visualization. The drugs were molecularly docked to the Hub genes, and the binding energy was calculated; the smaller the binding energy, the better and more stable the binding.

All bioinformatics statistical analyses were performed using R software, and P<0.05 was considered significant.

Experimental cells: Human umbilical vein endothelial cells (HUVECs) (Procell, Cat.CL-0675); Mouse RAW 264.7 cells (iCell, Cat.CL-0190).

Reagents: Oxidized low-density lipoprotein (ox-LDL) (Yiyuan Biotechnology, Cat.YB-002-1); Fetal bovine serum (Biosharp, Cat.BL205A); Trypsin-EDTA Solution (Biosharp, Cat.BL512A); RPMI-1640 medium (gibco, Cat.C11875500BT); DMEM medium (Pricella, Cat.PM150210); Penicillin/Streptomycin (Biosharp, Cat.BL505A); Oil red O (ORO) staining Kit (Solarbio, Cat.G1262); Paraformaldehyde, 4% (Solarbio, Cat.P1110); Cell count kit-8 (CCK-8) (Solarbio, Cat.CA1210); ROS detection Kit (BestBio, Cat. BB-47053); Mito-FerroGreen (DOJINDO, Cat.M489); Mounting medium, antifading (with DAPI) (Solarbio, Cat.S2110); Ferrostatin-1 (Fer-1) (APExBIO, Cat.A4371); JC-1 (Beyotime,Lot.C2003S); Serum-free cell freezing medium (Biosharp, Cat.BL203B); RNA extraction solution (Servicebio, Cat.G3013); Phosphate Buffered Saline (Servicebio, Cat.G4202); Chloroform substitute (Servicebio, Cat.G3014); RNA lysate (Servicebio, Cat.G3029); Water Nuclease-Free (Servicebio, Cat.G4700); SweScript All-in-One RT SuperMix for qPCR (One-Step gDNA Remover) (Servicebio, Cat.G3337); 2×Universal Blue SYBR Green qPCR Master Mix (Servicebio, Cat.G3326); isopropanol (Sinopharm, Cat.80109218); Anhydrous ethanol (Sinopharm, Cat.10009218).

Instruments: inverted fluorescence microscope (Nikon, Eclipse Ci); microplate reader (Tecan, Spark 10M); vortex mixer (Servicebio, SMV-3500); sealing instrument (Servicebio, FS-A20); microplate centrifuge (Servicebio, SMP-2); High-speed frozen microcentrifuge (DragonLab, D3024R); fluorescent quantitative PCR instrument (Bio-rad, CFX Connect); PCR instrument (Eastwin, ETC811).

Cells were cultured in RPMI-1640 medium (for HUVECs) or DMEM medium (for RAW 264.7 cells), supplemented with 10% fetal bovine serum and 1% penicillin/streptomycin and incubated at 37°C in humidified 5% CO₂. Cells were used from passages 3-8.

Cells were divided into three groups as follows: ①the Control group (CTRL): Cells were cultured normally without any treatment; ②the Model group (MOD): Cells were treated with ox-LDL (100 μg/mL) for 24 hours; ③the Fer-1 group (Fer-1): Cells were treated with ox-LDL (100 μg/mL) and Fer-1 (10 μM) for 24 hours.

Cells were fixed in 4% paraformaldehyde for 30 minutes and subsequently rinsed three times. Afterwards, the cells were stained with ORO for 1 hour at room temperature.

Cells were cultured in 96-well plates at a density of 2×10⁵ per well. After treatment with different interventions for 24 hours, 10 µL of CCK-8 was added to each well and incubated at 37°C for 1 hour. The absorbance was measured at 450 nm using a microplate reader.

DCFH-DA was diluted at 1:1000, the 12-well plates were removed from incubator, the medium was taken out. Then we washed cells with phosphate buffered saline for 3 times, 1 mL DCFH-DA was added, the cells were incubated at 37 °C for 2 hours, and then the cells were washed with serum-free medium for 3 times.

To detect intracellular mitochondrial Fe²⁺, Mito-FerroGreen were used according to the manufacturer’s protocol. The cells were inoculated in 12-well plates and cultured in a 5% CO₂ incubator at 37°C. After removing the medium, the cells were washed with serum-free medium for 3 times. Mito-Ferrogreen working solution 5μmol was added and cultured 1 hour in 5% CO₂ incubator at 37°C. After the supernatant was removed, the cells were washed with serum-free medium for 3 times.

JC-1 staining were used according to the manufacturer’s protocol. After the 12-well plates were removed from incubator, the culture medium was sucked out, the cells were washed once by phosphate buffered saline, and 0.5 mL cell culture medium and 0.5 mL JC-1 staining working solution were added, and fully mixed. The cells were incubated at 37 °C for 20 minutes in the incubator. After, the supernatant was removed and washed with JC-1 dyeing buffer twice. Add cell culture medium 1 mL.

Total RNA was extracted from cell lines using Trizol total RNA isolation reagent (Invitrogen) according to the manufacturer’s specifications and treated with Turbo DNase (Ambion). cDNA was synthesized from total RNA (0.5 mg) using random hexamers with the TaqMan cDNA Reverse Transcription Kit (Applied Biosystems). Primers were designed using Primer Express v3.0 software, and real-time PCR was performed using SYBR Select Master Mix (Applied Bio-Systems). All reactions were carried out on the 7500 Fast Real-Time PCR System (Applied Biosystem). The average of independent analyses for each gene and sample was calculated using the DD threshold cycle (Ct) method and was normalized to the endogenous reference control gene. The above primers were synthesized by Sangon Biotech (Shanghai) Co., Ltd (https://www.sangon.com/) ( Supplementary Table S2 ).

Fluorescence microscopy was used for observation and photography. Quantification was performed using Image J software, and results were presented as mean ± standard deviation. Statistical analysis was conducted with GraphPad Prism version 9.0 software. For data that met the criteria of normal distribution and homogeneity of variance, one-way analysis of variance (ANOVA) was used to compare differences among groups, followed by Tuke’s HSD post hoc test for pairwise comparisons; for data not meeting these criteria, the Kruskal-Wallis H non-parametric test was employed to compare differences between groups. A P<0.05 was considered statistically significant.

Differential gene expression analysis of the AS-related dataset GSE100927 identified 418 DEGs (adj.P<0.05 and |logFC|>1), comprising 295 up-regulated and 123 down-regulated genes ( Figures 2A–C ) ( Supplementary Table S3 ). To elucidate the functional implications of these DEGs in AS, we employed the GSEA and KEGG functional enrichment analysis method to uncover differential regulatory pathways between high- and low-expression groups. Our comprehensive analysis disclosed that the pathogenesis of AS chiefly encompasses organelle dysfunction (lysosomes, exosomes, endoplasmic reticulum, Golgi apparatus), cardiomyopathies (dilated cardiomyopathy, hypertrophic cardiomyopathy), metabolic dysregulation (tyrosine and propionic acid metabolism, valine, leucine, and isoleucine biosynthesis, glycosaminoglycan degradation), lipid metabolism disorder (lipids and AS, lipoprotein particle binding, cholesterol efflux, lipolytic metabolic processes), cell death mechanisms (ferroptosis, apoptosis), inflammatory processes (neutrophil extracellular trap formation, positive regulation of tumor necrosis factor, IL-6, IL-8, IL-12, and IL-1β production), and immune dysregulation (primary immunodeficiency, natural killer cell-mediated cytotoxicity, B-cell receptor signaling pathways). These pathways are closely associated with various organelles, such as lysosomes, endoplasmic reticulum (ER), and Golgi apparatus ( Figures 2D–H ).

From the FerrDb database, we retrieved 484 FeRGs ( Supplementary Table S4 ), while 5151 LMRGs were obtained from public databases and literature queries ( Supplementary Table S5 ). The intersection of these datasets with DEGs yielded six ASLMRFeGs (CTSB, CYBB, DPP4, HILPDA, HMOX1, IL1B) ( Figure 3A ). Among these, CTSB (logFC=1.530), CYBB (logFC=1.116), DPP4 (logFC=1.430), HMOX1 (logFC=1.839), and IL1B (logFC=1.447) exhibited higher expression in AS, whereas HILPDA (logFC=-1.172) displayed lower expression ( Figures 3B, C ). In addition, CTSB, CYBB, DPP4, HILPDA, and IL1B are ferroptosis driver genes, whereas HMOX1 is a bidirectional regulator of ferroptosis (driver and suppressor). To delve deeper into the roles of ASLMRFeGs in AS progression, we conducted correlation analyses. The results illuminated robust synergistic effects among CTSB, CYBB, DPP4, HMOX1, and IL1B, whereas HILPDA exhibited an antagonistic relationship with CTSB, CYBB, DPP4, and HMOX1 ( Figure 3D ).

Expanding our investigation into the immune microenvironment of AS, we analyzed specific immune cell types infiltrating AS tissues using the CIBERSORT method ( Figure 3E ). Our findings revealed significantly elevated levels of memory B cells, activated memory CD4 T cells, follicular helper T cells, γ-δ T cells, M0 macrophages, and activated mast cells in AS. Conversely, significantly lower levels of naive B cells, plasma cells, CD8 T cells, resting memory CD4 T cells, regulatory T cells, resting NK cells, monocytes, and resting mast cells were observed in AS ( Figure 3F ). Furthermore, correlation analysis unveiled significant positive associations of CTSB, CYBB, DPP4, and HMOX1 with M0 macrophages, along with significant negative correlations with M1 macrophages, resting mast cells, and resting memory CD4+ T cells. IL1B showed a negative association with M2 macrophages and resting mast cells but a positive link with activated mast cells ( Figure 3G ). These findings suggest that ASLMRFeGs contribute to the regulation of the immune microenvironment in AS by influencing specific immune cell infiltration and activation states.

To classify the 69 AS samples based on the expression profiles of the six ASLMRFeGs, a consensus clustering algorithm was employed, resulting in the identification of two distinct and stable groups: the C1 cluster (n=33) and the C2 cluster (n=36) ( Figures 4A–E ). Further analysis of the expression differences of the six ASLMRFeGs between the C1 and C2 clusters revealed that CTSB, CYBB, DPP4, HMOX1, and IL1B were highly expressed in the C2 cluster, while HILPDA was highly expressed in the C1 cluster ( Figures 4F, G ).

Immune infiltration analyses of the two clusters demonstrated that the C1 cluster had a significantly higher abundance of memory B cells, naive CD4 T cells, resting memory CD4 T cells, M1 macrophages, and resting mast cells. In contrast, the C2 cluster exhibited significantly higher levels of follicular helper T cells and M0 macrophages ( Figures 4H, I ).

GSVA was then utilized to investigate potential biological and functional differences between the two clusters. The results indicated that the C2 cluster was primarily enriched in pathways related to inflammation and immune response, including leukocyte proliferation, regulation of dendritic cell chemotaxis, neutrophil chemotaxis, antigen processing and presentation, the Toll-like receptor signaling pathway, and the NOD-like receptor signaling pathway. Conversely, the C1 cluster was mainly associated with pathways affecting the structure and function of the cardiovascular system and cell signaling, such as glycoprotein complex, regulation of cell communication by electrical coupling involved in cardiac conduction, dilated cardiomyopathy, hypertrophic cardiomyopathy, arrhythmogenic right ventricular cardiomyopathy, vascular smooth muscle contraction, the TGF-beta signaling pathway, and the Wnt signaling pathway ( Figures 4J, K ). These findings suggest that, based on the expression of ASLMRFeGs, AS samples can be divided into two subgroups with significantly different biological functions, particularly in terms of inflammation and immune response.

By utilizing WGCNA, co-expression networks and modules were constructed for both control and AS samples ( Figures 5A, B ). Upon setting the soft threshold to 6, 8 distinct modules with different colors were identified, among which the “turquoise” module exhibited the strongest positive correlation with AS (r=0.72) ( Figures 5C-E ). This particular module comprised 6069 module genes and 1425 module Hub genes ( Figure 5F ). Moreover, WGCNA analysis revealed module Hub genes that were associated with lipid metabolism pathway ferroptosis in AS samples ( Figures 5G, H ). When the soft threshold was set to 11, a total of 8 modules with distinct colors were detected. Among them, the “black,” “blue,” “brown,” “green,” and “turquoise” modules are significantly correlated with the C2 cluster, and once again, the “turquoise” module showcased the highest positive correlation with the C2 cluster (r=0.87) ( Figures 5I–K ). In order to screen the most relevant molecular modules for lipid metabolism and ferroptosis, we enriched the genes in the different modules and showed that the “turquoise” module was closely related to inflammatory response, immune response, lipid metabolism, and ferroptosis; the “black” module was related to cell cycle; the “blue” module was related to muscle contraction, regulation of actin cytoskeleton, and cell signaling; and the “green” module was related to cell adhesion and angiogenesis, while the “brown” module was not significantly enriched. Therefore, the “turquoise” module, which includes 3284 module genes and 1515 module Hub genes, was used for the study ( Figure 5L ) ( Supplementary Table S6 ).

The intersection of module genes from the two “turquoise” modules and DEGs resulted in 225 candidate module Hub genes ( Figure 6A ), and constructed a PPI network with 118 nodes and 359 edges ( Figure 6B ). The analysis of these genes using gene ontology and KEGG revealed their involvement in various biological functions and pathways. The gene ontology analysis encompassed three main categories: biological process, cellular component, and molecular function. Among the 225 candidate module Hub genes, 189 were linked to biological processes, 55 to cellular components, and 41 to molecular functions, while the KEGG enrichment analysis identified 52 signaling pathways (P<0.05). The biological process analysis indicated the participation of these genes in several processes, including inflammatory response, immune response, signal transduction, the positive regulation of tumor necrosis factor, macrophage cytokine, IL-6, IL-8, and IL-10 production, and lipoprotein catabolic process. The cellular component analysis highlighted the role of various cellular membranes and organelles, especially the lysosomal membrane, late endosome membrane, endocytic vesicle membrane, endosome membrane, ER-to-Golgi transport vesicle membrane, Golgi membrane, and mitochondrial outer membrane. Furthermore, the molecular function analysis emphasized their involvement in MHC class II protein complex binding, transmembrane signaling receptor activity, superoxide-generating NAD(P)H oxidase activity, signaling receptor activity, actin filament binding, and cysteine-type endopeptidase activator activity involved in apoptotic processes ( Figures 6E, F ). The pathway analysis revealed enrichment in various pathways, including the B cell receptor signaling pathway, the Fc epsilon RI signaling pathway, the chemokine signaling pathway, the NOD-like receptor signaling pathway, the Toll-like receptor signaling pathway, the T cell receptor signaling pathway, the Rap1 signaling pathway, Th17 cell differentiation, natural killer cell-mediated cytotoxicity, complement and coagulation cascades, lysosome, phagosome, efferocytosis, neutrophil extracellular trap formation, Fc gamma R-mediated phagocytosis, leukocyte transendothelial migration, and lipid and AS ( Figure 6G ). The above results show that the pathogenesis of ferroptosis in the AS lipid metabolism pathway involves abnormalities in several biological processes, including cellular communication, immune response, and metabolic response.

To further screen the candidate Hub genes, 225 candidate module Hub genes were externally validated by single-gene differential analysis in the GSE43292, GSE28829, and GSE57691 datasets. This validation revealed that 76 module Hub genes were significantly differentially expressed across all four datasets ( Figures 6C, D ) ( Supplementary Table S7 ). These 76 genes were then imported into the DAVID database for gene ontology and KEGG pathway analysis. The results showed that 61 genes were correlated with biological processes, 36 with cellular components, and 14 with molecular functions. Additionally, KEGG enrichment analysis identified 38 signaling pathways (P<0.05).

The biological process analysis indicated involvement in several processes, including immune response, inflammatory response, signal transduction, positive regulation of phagocytosis, actin filament organization, positive regulation of tumor necrosis factor and interleukin-6 production, superoxide anion generation, leukocyte cell-cell adhesion, negative regulation of neuron apoptotic process, positive regulation of macrophage fusion, and interleukin-10-mediated signaling pathway. The cellular component analysis identified the phagolysosome, late endosome membrane, Golgi membrane, endosome membrane, ER-to-Golgi transport vesicle membrane, lysosome, phagocytic vesicle, and lamellipodium. The molecular function analysis highlighted amyloid-beta binding, MHC class II protein complex binding, actin binding, signaling receptor activity, cytokine binding, and superoxide-generating NADPH oxidase activator activity ( Figures 6H, I ).

Pathway analysis indicated enrichment in the Rap1 signaling pathway, the Fc epsilon RI signaling pathway, phagosome, neutrophil extracellular trap formation, leukocyte transendothelial migration, cell adhesion molecules, Th1 and Th2 cell differentiation, Th17 cell differentiation, lipid and AS, and natural killer cell-mediated cytotoxicity. These findings suggest that AS may modulate lipid metabolism and ferroptosis by regulating inflammation and immune responses ( Figure 6J ).

Subsequently, these 76 module Hub genes were imported into the STRING database to construct a PPI network with a “score” ≥ 0.7, which had 40 nodes and 79 edges ( Figure 6D ). 8 candidate Hub genes (TYROBP, CSF1R, PTPRC, ITGAM, LCP2, CYBB, HLA-DRA, and C1QA) with degree centrality, betweenness centrality, and closeness centrality all greater than the mean were filtered out ( Supplementary Table S8 ), and the constructed PPI network of candidate Hub genes includes 8 nodes and 15 edges ( Figure 7A ). To further refine the selection of characterized genes, 3 types of machine learning algorithms were employed: LASSO regression, SVM-RFE, and RF. LASSO regression identified 6 key genes (TYROBP, CSF1R, PTPRC, LCP2, HLA-DRA, and C1QA) ( Figures 7B, C ), SVM-RFE identified 6 featured genes (CSF1R, PTPRC, TYROBP, LCP2, C1QA, and ITGAM) ( Figure 7D ), and RF obtained 4 important genes (CSF1R, TYROBP, C1QA, and LCP2) ( Figures 7E, F ) ( Supplementary Table S9 ). Ultimately, an intersection of the genes identified by these methods yielded 4 characterized genes: TYROBP, CSF1R, LCP2, and C1QA, which were all up-regulated in AS ( Figure 7G ).

Similarly, external validation of single-gene differential analysis of six ASLMRFeGs in the GSE43292, GSE28829, and GSE57691 datasets showed that CYBB, HMOX1, and IL1B were significantly differentially expressed in these datasets and up-regulated in AS. Therefore, we considered these three ASLMRFeGs and four characterized genes as candidate Hub genes, namely, CYBB, HMOX1, IL1B, CSF1R, TYROBP, C1QA, and LCP2 ( Figures 7H–J ).

We validated this using the RAW 264.7 and HUVECs dual cell lines, one belonging to endothelial cells and one to macrophages. Firstly, we needed to confirm that ox-LDL-induced HUVECs and RAW 264.7 cells clearly caused abnormal lipid metabolism, so ORO staining was used to detect intracellular lipid accumulation, and CCK-8 was used to detect cell activity. The results showed that there was only a small amount of ORO staining observed in normal HUVECs and RAW 264.7 cells, and extensive ORO staining was visible in the cytoplasm of ox-LDL-treated HUVECs and RAW 264.7 cells in the MOD group, with increased cell volume and rounded or irregular morphology. The relative area of ORO was significantly increased (P<0.01) and cell activity was significantly decreased (P<0.01) in the MOD group compared with the CTRL group, suggesting that ox-LDL-induced HUVECs and RAW 264.7 cells resulted in an intracellular abnormality of lipid metabolism, which led to significant lipid accumulation and decreased cell activity in the cells ( Figures 8A–F ).

To verify whether abnormal lipid metabolism could promote ferroptosis, we used Mito-FerroGreen to detect mitochondrial Fe²⁺ content. The results showed that, compared with the CTRL group, the MOD group significantly elevated mitochondrial Fe²⁺ content (P<0.01), suggesting that ox-LDL induces ferroptosis in HUVECs and RAW 264.7 cells and that abnormalities in lipid metabolism promote ferroptosis ( Figures 8G–J ). Then, we used Fer-1 (an inhibitor of ferroptosis) for subsequent experiments, and the results showed that the Fer-1 group significantly alleviated intracellular lipid accumulation (P<0.01), elevated cell activity (P<0.05), and reduced mitochondrial Fe²⁺ content (P<0.01) compared with the MOD group, suggesting that inhibition of ferroptosis improves lipid metabolism and cell activity to a certain extent ( Figures 8A–F ).

In order to predict the potential mechanism of lipid metabolism in AS and ferroptosis, we used qRT-PCR to examine the expression levels of ASLMRFeGs and Hub genes. The results showed that ox-LDL treatment (MOD group) significantly up-regulated the expression of CYBB, HMOX1, IL1B, CSF1R, TYROBP, and LCP2 genes, while down-regulating C1QA gene compared with the CTRL group in both HUVECs and RAW 264.7 cells. Treatment with the ferroptosis inhibitor Ferrostatin-1 (Fer-1 group) significantly attenuated this ox-LDL-induced upregulation of CYBB, HMOX1, IL1B, CSF1R, and TYROBP genes compared with the MOD group in both cell lines. Additionally, Fer-1 attenuated LCP2 upregulation specifically in RAW 264.7 cells ( Figures 8K, L ).

The previous bioinformatics analysis results indicated that the mechanisms by which lipid metabolism regulates AS and ferroptosis are closely related to various organelles, particularly mitochondria, lysosomes, the ER, and the Golgi apparatus. Therefore, we investigated the role of mitochondria using three different methods: the JC-1 fluorescent dye to measure the relative fluorescence intensity ratio (red/green fluorescence) of mitochondrial membrane potential, with a higher ratio indicating better mitochondrial function; DCFH-DA fluorescent staining to assess the content of ROS; and qRT-PCR to detect the relative expression level of mtDNA-ND1 in mitochondria to determine mitochondrial function.

The results showed that the fluorescence intensity ratio of mitochondrial membrane potential was significantly weaker, the fluorescence intensity of ROS was significantly stronger, and the relative expression level of ND1 in mitochondria was decreased in the MOD group compared with the CTRL group. This suggests that ox-LDL-induced mitochondrial dysfunction in HUVECs ( Figures 9A–E ) and RAW 264.7 cells ( Figures 9F–J ) promote ROS production. However, Fer-1 antagonized the mitochondrial-damaging effects induced by ox-LDL, indicating that inhibition of ferroptosis improved mitochondrial dysfunction induced by abnormalities in lipid metabolism.

Next, we found that 6 hours of ox-LDL intervention upregulated the relative expression of LAMP1 and LAMP2 genes in HUVECs and the LAMP2 gene in RAW 264.7 cells. At 12 hours, ox-LDL intervention continued to upregulate LAMP2 gene expression in HUVECs. However, after 24 hours of ox-LDL intervention, both LAMP1 and LAMP2 gene expression were downregulated in both HUVECs ( Figures 9K, L ) and RAW 264.7 cells ( Figures 9N, O ). In contrast, Fer-1 treatment increased the expression of both LAMP1 and LAMP2 genes. These findings suggest that abnormal lipid metabolism may induce lysosomal activation in the early stages, followed by lysosomal damage in later stages. Inhibition of ferroptosis may help mitigate lysosomal damage.

GOLGA2 is a classical marker of the Golgi apparatus. We found that the MOD group downregulated GOLGA2 gene expression compared with the CTRL group, whereas Fer-1 upregulated GOLGA2 gene expression compared with the MOD group ( Figures 9M, P ).

The ER is an important organelle for lipid and protein synthesis. Under various stimuli, overloaded protein synthesis exceeds the storage capacity of the ER, promoting ERS. GRP78 and CHOP are commonly used as markers of ERS. We also detected genes of the three pathways of ERS (IRE1, ATF6, PERK, XBP1, and EIF2A). The results indicated that in HUVECs, the MOD group significantly increased the expression of ATF6, GRP78, CHOP, and EIF2A genes while decreasing the expression of IRE1, PERK, and XBP1 genes compared to the CTRL group. In contrast, the Fer-1 group upregulated IRE1, ATF6, PERK, and XBP1 genes and downregulated GRP78 and CHOP gene expression compared to the MOD group. In RAW 264.7 cells, the MOD group upregulated IRE1, PERK, GRP78, CHOP, and EIF2A genes while downregulating ATF6 and XBP1 genes compared to the CTRL group. In contrast, the Fer-1 group downregulated IRE1, PERK, GRP78, CHOP, and EIF2A genes compared to the MOD group ( Figures 9Q, R ). This suggests that abnormal lipid metabolism significantly activates ERS, and inhibition of ferroptosis suppresses ERS.

Our experiments demonstrated that abnormal lipid metabolism can damage mitochondrial mtDNA, causing mitochondrial dysfunction, promoting ROS production, and increasing mitochondrial Fe²⁺ content, which in turn promotes ferroptosis. Mitochondria, as a major source of ROS in cells, play a crucial role in regulating cellular metabolism, signal transduction, and death signaling, especially in iron metabolism as well as material and energy metabolism (53). Mitochondrial redox function plays a decisive role in the development of AS. By enhancing mitochondrial oxidative metabolism promotes fatty acid degradation, reduces intracellular lipid accumulation, limits foam cell formation, and delays the onset of AS (54, 55), whereas inhibition of ferroptosis can attenuate AS by decreasing lipid peroxidation and endothelial dysfunction (56). Studies have revealed distinct changes in mitochondria during ferroptosis, including marked shrinkage, increased membrane density, and a reduction or disappearance of mitochondrial cristae, which is not consistent with characteristics of other cell death, such as apoptosis and cell necrosis (57). To counteract the detrimental effects of ferroptosis on mitochondria, mitochondria significantly inhibit oxidative stress-induced ferroptosis by overexpressing mitochondrial ferritin or mitochondrial catalase, which inhibit mitochondrial iron overload or ROS accumulation (58). In addition, mitochondria play a central role in fatty acid metabolism by providing specific lipid precursors for lipid oxidation, which is a key step in cellular ferroptosis (53). Mitochondrial energy metabolism is also closely intertwined with ferroptosis, and inhibition of the mitochondrial tricarboxylic acid cycle or the electron transfer chain can attenuate mitochondrial membrane potential hyperpolarization, lipid peroxide accumulation, and ferroptosis (59).

Lysosomes, as acidic membrane-bound organelles, serve as the endpoints of various vesicular transport pathways, including endocytosis, phagocytosis, and autophagy pathways, which are essential for promoting lipid catabolism and transport and maintaining cellular lipid homeostasis (60). Lysosomal proteins LAMP1 and LAMP2 have been shown to bind cholesterol with high affinity and specificity through their luminal domain 1 and facilitate the flow of cholesterol through the NPC2-NPC1 export pathway (61). Additionally, lysosomes may rely on membrane contacts with other organelles for lipid transfer, and cholesterol transfer in the ER-to-lysosomal direction was enhanced in a VAP-dependent manner by overexpression of STARD-3 (62); STARD-3 also transfers cholesterol from lysosomes to mitochondria to provide precursors for steroid hormones synthesized in mitochondria (63).

Recently, it was found that inhibition of lysosome-dependent cell death limited Erastin-induced ferroptosis, suggesting that ferroptosis is a lysosomal cell death process and that lysosomes contribute to ferroptosis through mechanisms such as autophagy activation, the release of lysosomal cathepsins, and the accumulation of lysosomal iron or nitric oxide (64, 65). Many studies have shown that excessive activation of autophagy and lysosomes degrades ferritin, increases unstable iron accumulation, and promotes iron-dependent cell death (66). Cathepsin B (CTSB) is a mediator of organelle-specific initiation of ferroptosis, and activation of the transcription factor STAT3 upregulates CTSB, which is translocated from the lysosome to the nucleus, leading to DNA damage and subsequent STING1-dependent ferroptosis (65). Conversely, CTSB-dependent albumin catabolism promotes glutathione synthesis by exporting cystine from lysosomes via the transporter cystinosis fuels, thereby inhibiting lethal lipid peroxidation (67). Additionally, it was found that lysosomal damage caused by mitochondrial dysfunction induced ferroptosis (68). This indicates that lysosomes regulate lipid metabolism through membrane contacts with other organelles, influencing autophagy, cathepsin release, and thus ferroptosis. Our study suggests that endothelial cell injury and foam cell models may lead to activation of lysosomes and promotion of autophagy in the initial phase, followed by lysosomal over-excitation, causing dysfunction, whereas inhibition of ferroptosis improves lysosomal function.

The ER plays a critical role in the synthesis, processing, packaging, and transportation of proteins and lipids. When unfolded or misfolded proteins accumulate in the ER, it leads to a pathological condition known as ERS. To counter ERS, cells initiate the unfolded protein response (UPR) signaling cascade by activating the expression of molecular chaperones, regulating lipid synthesis, and promoting ER-associated degradation, which helps reduce the levels of unfolded or misfolded proteins and restore ER homeostasis (69, 70). ERS has been identified as a key regulator of cholesterol deposition, macrophage differentiation, and endothelial dysfunction. The ER is directly involved in the formation of lipid droplets and the maintenance of lipid homeostasis. When disorders of lipid metabolism are exacerbated, it promotes ERS and URP, disrupting the normal metabolism of adipose tissue, triggering an inflammatory response in adipocytes, enhancing the secretion of adipokines, and causing ectopic lipid deposition, which contributes to the development of metabolic diseases (71). This situation is particularly pronounced in obese patients, where overnutrition stimulates adipogenesis, creating an imbalance in the body’s microenvironment, inducing progressive hypertrophy of adipocytes, stimulating the ER to synthesize more proteins for the formation of lipid droplets, and altering the flow and thickness of ER membranes via the ER sensors IRE1α and PERK, which lead directly to ERS (72). Furthermore, ERS is an important cause of endothelial and macrophage apoptosis in advanced lesions (73). Inhibiting macrophage ERS promotes the polarization of differentiated M2 macrophages toward an M1 phenotype, thereby inhibiting foam cell formation (74). ERS also contributes to an imbalance between nitric oxide and ROS, causing oxidative stress, damaging endothelial cells, and promoting the progression of AS and plaque formation (73, 75).

The classical UPR signaling pathway involves the activation of three classes of ER membrane proteins: IRE1, PERK, and ATF6 (71). In normal conditions, the chaperone protein Bip and GRP78 bind to these three proteins, preventing their activation and keeping them in an inactive state. However, when external stimuli or changes in ER homeostasis result in the excessive accumulation of misfolded or unfolded proteins, Bip binds to the abnormal proteins in the ER lumen, dissociates from and activates IRE1α, PERK, and ATF6 to carry out their respective functions. This activation enhances ER protein folding capacity, inhibits intracellular protein synthesis, and restores ER homeostasis (76). Nevertheless, if the degree of ERS surpasses the UPR’s regulatory range, it disrupts ER homeostasis, leading to various metabolic abnormalities associated with ER function. GRP78 and CHOP are commonly used as markers of ERS, and our experiments also validated the three branches of ERS. First, based on the expression of GRP78 and CHOP genes, it was shown that abnormal lipid metabolism significantly activated ERS, which might be excessive, and the use of Fer-1 inhibited the excessive ERS of the cells. It was further found that the abnormal lipid metabolism regulated the UPR branch of IRE1 and PERK, and the use of Fer-1 regulated the branches of IRE1, PERK, and ATF6.

In the first branch, autophosphorylated IRE1 activates the RNase structural domain and catalyzes the splicing of XBP1 mRNA, generating the active transcription factor XBP-1s. XBP-1s is involved in inflammation, cell survival, lipid metabolism, and calcification to regulate endothelial cell proliferation, transformation, and apoptosis, smooth muscle cell phenotypic switching, and the accumulation of foam cells that thereby affect AS (77). Furthermore, the overexpression of IREI and XBP1 has been shown to increase cellular susceptibility to ferroptosis (78, 79). Additionally, IRE1 activates c-Jun NH2 terminal kinase (JNK) in response to ERS. Inhibition of the IRE1/JNK pathway not only reduces ERS-induced apoptosis and improves vascular endothelial dysfunction (80) but also attenuates ferroptosis in acute kidney injury (81). Collectively, these studies indicate the IRE1 branch plays a crucial role in AS progression through XBP-1s/JNK signaling, ERS-induced apoptosis, and enhanced ferroptosis susceptibility. Unexpectedly, our experiments revealed that ox-LDL-induced lipid accumulation promotes ERS while concomitantly downregulating XBP1 expression. This apparent paradox, ERS induction with reduced IRE1/XBP1 signaling, suggests a potential compensatory adaptation may occur under sustained lipid stress. We hypothesize that persistent lipid overload triggers feedback inhibition of the IRE1/XBP1 pathway to mitigate excessive cellular stress responses, which could consequently reduce ferroptosis susceptibility. This interpretation aligns with our observation that ferroptosis inhibition (Fer-1) alleviates both lipid accumulation and ERS severity. Under such conditions of reduced lipid stress, restoring IRE1/XBP1 signaling could theoretically help reestablish ER homeostasis through ER-associated degradation of misfolded proteins and upregulation of pro-homeostatic genes (77), potentially delaying AS progression. However, the precise mechanisms underlying lipid-mediated suppression of IRE1/XBP1 and its functional consequences require further experimental investigation.

In the second branch, PERK plays a dual role, serving not only to reduce the protein load of ERS but also to signal the cell death pathway. Phosphorylated PERK phosphorylates eIF2α, thereby inhibiting mRNA translation, reducing protein influx into the ER, and facilitating the translation of ATF4. CHOP is an important target gene regulated by ATF4, and sustained ERS causes CHOP to activate genes responsible for encoding apoptosis (82). In addition, PERK acts as a negative regulator of ferroptosis, and the occurrence of ferroptosis is accompanied by the activation of the PERK-eIF2α-ATF4-CHOP signaling pathway, and inhibition of ERS reduces ferroptosis (83–85). Our experiments found that abnormal lipid metabolism up-regulated the expression of the EIF2A and CHOP genes; the use of Fer-1 improved the expression of the PERK, EIF2A, and CHOP genes, possibly due to the negative feedback regulation of up-regulated EIF2A to inhibit the expression of PERK.

In the third branch, ATF6 is released and translocated to the Golgi, where it undergoes cleavage by resident proteases (site 1 and site 2 proteases) and subsequently moves to the nucleus to activate the transcription of various chaperone molecules. Previous studies have shown that increased expression of ATF6α contributes to enhanced resistance against ferroptosis (86).

Additionally, the Golgi apparatus plays a key role in regulating protein and lipid synthesis, modification, and distribution within the cell. Golgi stress not only triggers the production of ROS but also promotes the accumulation of lipid peroxides, ultimately inducing ferroptosis (87). In short, the module genes obtained by WGCNA provide a preliminary indication of the mechanisms by which abnormalities of lipid metabolism regulate AS and ferroptosis, in which inflammation and immune responses play key roles. In particular, the abnormalities of lipid metabolism may affect ferroptosis and AS by regulating cell organelle function.