AS-related expression profile data were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The GSE43292 dataset (arterial tissues from 32 controls and 32 AS patients) was used as the training set and the GSE100927 dataset (arterial tissues from 12 controls and 29 AS patients) was used as the validation set. In total, 259 ferroptosis-related genes (Ferr-RGs) were obtained from the FerrDb database (http://www.zhounan.org/ferrdb/current/) (16) (Supplementary Table S1) and 24 m6A regulators were obtained from the literature (17) (Supplementary Table S2).

To identify the differentially expressed genes (DEGs) between different samples, the data in GSE43292 were analyzed for DEGs associated with AS using the R language “limma” package (version 3.50.1) (18) with an adjusted p-value (p_adj) <0.01 as a screening criterion. Volcano and heat maps were plotted using “ggplot2” (version 3.3.5) and “pheatmap” (version 1.0.12), respectively.

The intersection of DEGs with m6A-related genes (m6A-RGs) was taken to obtain the differentially expressed m6A-associated genes (DE-m6A-RGs), using Venn (version 1.11) (19) to create the Venn diagrams. Similarly, the differentially expressed ferroptosis-related genes (DE-Ferr-RGs) were obtained. Pearson’s correlation was calculated between these genes. Genes that met the criteria (p < 0.001 and |R|>0.4) were recognized as differentially expressed m6A-associated ferroptosis-related genes (DE-m6A-Ferr-RGs), and the results of the calculation were visualized in a heat map.

To identify genes associated with different traits, WGCNA (version 1.70.3) (20) was performed on the samples from GSE43292. A hierarchical clustering tree was first constructed for all the samples (N = 64), and a sample dendrogram and a trait (AS vs. control) heat map were constructed. Next, the appropriate soft threshold was selected based on the near-scale-free topological criteria to construct gene modules. The modules were then gathered according to the criteria of the dynamic tree-cut algorithm. The correlations of each gene module with traits (control and AS) were calculated, the conditions for selecting key modules were: non-gray modules, p-values <0.05, and |cor| >0.5. Finally, module membership (MM) and gene significance (GS) screening were performed, and the key module genes were screened using the criteria of |GS| >0.4 and |MM| >0.8, respectively.

The DE-m6A-Ferr-RGs were intersected with the key module genes to obtain candidate m6A-Ferr-related signature genes. To investigate the biological pathways involving the candidate m6A-Ferr-related signature genes, the Gene Ontology (GO) functional enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses of the candidate m6A-Ferr-related signature genes were performed using the “clusterprofiler” package (version 4.2.2) (21).

To further identify m6A-Ferr-related signature genes, LASSO analysis (22) of the candidate m6A-Ferr-related signature genes was performed using glmnet (version 4.1.2), with the “family” parameter set as “binomial” and the cross-validation parameter “nfolds” adjusted to 10 to obtain the cross-validated error maps and gene coefficient maps. Thus, the m6A-Ferr-related signature genes were finally obtained. Subsequently, expression analysis of the m6A-Ferr-related signature genes was implemented in the GSE43292 and GSE100927 datasets, respectively. Receiver operating characteristic (ROC) curves were drawn using the pROC package (version 2.3.0) (23) to evaluate the diagnostic value of the m6A-Ferr-related signature genes. The larger the area under the curve (AUC), the higher the accuracy. The R package “rms” (version 6.2-0) (24) was used to construct the nomogram to predict the incidence of AS. Calibration curves were plotted using “regplot” (version 1.1) (25). Decision curves (DCA) was also plotted using the “rmda” package.

To find the significant pathways for the m6A-Ferr-related signature genes, single-gene gene set enrichment analysis (GSEA) (26) was performed. The filtering criteria were |NES| >1, p_adj <0.05, and q-value <0.25.

To investigate whether there were differences in immune cells between the atherosclerosis and control groups, we used the false discovery rate (FDR) correction to calculate the proportion of 22 kinds of immune cells in the AS samples. Heat maps and boxplots were drawn for visualization. Lollipop mapping of the 22 immune cells and m6A-Ferr-related signature genes were plotted to show the correlation analysis result.

The relative RNA expression level genes were quantified by a real-time polymerase chain reaction procedure. Thus, 12 frozen tissue samples were obtained from the College of Forensic Medicine, Kunming Medical University, of which samples 1–5 were the control group and 6–12 were the AS group. This study was approved by the Medical Ethics Committee of the Second Affiliated Hospital of Kunming Medical University in China. All patients had signed an informed consent form. The expression of m6A-Ferr-related signature genes was further validated via RT-qPCR. The total RNA of the 12 samples was extracted using TRIzol (Ambion, Austin, USA) according to the manufacturer's instructions. Reverse transcription of the total RNA to cDNA was carried out using a SureScript-First-strand-cDNA-synthesis-kit (Servicebio, Wuhan, China) based on the manufacturer's instructions. RT-qPCR was performed utilizing the 2xUniversal Blue SYBR Green qPCR Master Mix (Servicebio, Wuhan, China). The primer sequences for the PCR are shown in Supplementary Table S3. GAPDH was used as an internal reference gene. The 2^−ΔΔCt method was utilized to calculate the expression of key genes (27).

All bioinformatics analyses were undertaken in R language. The rank sum test was employed to contrast the data from different groups. p < 0.05 was considered to represent a significant difference.

The differential analysis obtained a total of 6,167 DEGs in the AS and control samples from the GSE43292 dataset, of which 3,112 genes were upregulated and 3,055 genes were downregulated in the AS samples (Figure 1A). The top 100 upregulated and downregulated genes are presented in a heat map (Figure 1B). The intersection of the DEGs with 24 m6A-RGs and 259 Ferr-RGs resulted in nine DE-m6A-RGs and 104 DE-Ferr-RGs, respectively (Figure 1C). The results of the correlation between the DE-m6A-RGs and the DE-Ferr-RGs genes are shown in Figure 1D. In total, 113 genes satisfied p <0.001 and |R| >0.4 and were considered as DE-m6A-Ferr-RGs for subsequent analysis.

WGCNA was performed in the GSE43292 dataset. The hierarchical clustering tree results, as shown in Figure 2A, showed no significant outlier samples, so all samples were used for the subsequent analysis. The clustering and trait heat map of the AS and control samples is shown in Figure 2B. As shown in Figure 2C, nine with the best soft thresholds were chosen to construct the gene modules, and 23 gene modules were obtained, among which the gray module was non-sense (Figure 2D). The heat map shows the correlation between the 23 gene modules and the AS and control groups (Figure 2E). The key modules associated with AS were the blue, tan, green, yellow, brown, and green modules according to a p-value of <0.05 and |cor| of >0.5 (Figure 2F). The key module genes were screened using the criteria of |GS| >0.4 and |MM| >0.8, respectively, and a total of 2,373 module genes were obtained (Figure 2G). Finally, the DE-m6A-Ferr-RGs were intersected with the key module genes to obtain 48 candidate m6A-Ferr-related signature genes (Figure 2H).

To determine the potential biological roles of the selected candidate m6A-Ferr-related signature genes, we performed enrichment analyses. The GO functional enrichment analysis resulted in a total of 295 items, including 250 biological process (BP) items (cellular response to chemical stress, response to oxidative stress, neutrophil degranulation, etc.), 31 cellular components items (tertiary granule, endocytic vesicles, membrane rafts, etc.), and 14 molecular functions items (ubiquitin protein ligase binding, ferrous iron binding, dioxygenase activity, etc.) (Supplementary Table S4), and Figure 3A shows the top 10 ranked items under each GO classification. The results of the KEGG pathway analysis showed that a total of 11 pathways were enriched, and the enriched pathways [ferroptosis, autophagy-animal, chemical carcinogenesis-reactive oxygen species (ROS), etc.] are shown in Figure 3B.

Screening for m6A-Ferr-related signature genes was conducted using LASSO regression analysis. The LASSO analysis screened six genes (AGPAT3, NOX4, CDO1, CYBB, ATG7, and SLC3A2) that were recorded as m6A-Ferr-related signature genes (Figures 4A,B). ROC curves of the LASSO model with AUC values of 0.880 and 0.741 in the GSE43292 dataset and the GSE100927 dataset, respectively (Figures 4C,D), indicated that the model performed well. Thus, a nomogram was plotted based on the six genes. The adjusted C-index of the nomogram was 0.835, indicating that the selection of the genes was appropriate (Figure 4E). The calibration and decision curves confirmed the above conclusion (Figures 4F,G).

The expression of the m6A-Ferr-related signature genes in the AS and control samples from GSE43292 is shown in Figure 5A, in which AGPAT3, ATG7, CYBB, and SLC3A2 were highly expressed and NOX4 and CDO1 were slightly expressed in AS group. Furthermore, the expression of the m6A-Ferr-related signature genes was verified in GSE100927. As shown in Figure 5C, the expression trend of the genes in the GSE100927 dataset was completely consistent with the GSE43292 dataset. In addition, in the GSE43292 and GSE100927 datasets, the AUC values of all six m6A-Ferr-related signature genes were greater than 0.7, indicating that they could distinguish AS patients from the control samples and had potential diagnostic value (Figures 5B,D). We collected 12 histological samples for RT-qPCR validation, and the detailed disease information of the corresponding patient samples is presented in Supplementary Table S5. The RT-qPCR results revealed that the expression of SLC3A2, NOX4, and CDO1 was consistent with the transcriptome level. The expression trend of AGPAT3, CYBB, and ATG7 was consistent but there was no difference between the AS and control groups (Figure 5E).

To understand the biological function and the involved signaling pathways of the m6A-Ferr-related signature genes, GSEA was performed. The top 10 GO entries and KEGG pathways of the six genes are displayed in Figures 6A–F, 7A–F. We found that these genes were involved in B-cell-mediated immunity, activated immune response, adaptive immune response, the T-cell receptor signaling pathway, the B-cell receptor signaling pathway, and autoimmune thyroid disease among other BP and pathways (Figures 6A–F, 7A–F).

Since the GSEA results suggested that the m6A-Ferr-related signature genes in AS were linked to immune-related functions, we then performed an immune infiltration analysis. By analyzing the ratio of 22 immune cells in the AS and control samples, an abundance map of immune cell infiltration proportions was constructed to demonstrate the immune cell content in the samples (Figure 8A). A correlation heat map and p-value correlation heat map of the immune cell ratio are shown in Figures 8B,C. A total of two types of immune cells (CD8+ T cells and activated CD4+ memory T cells) showed significant differences between the AS and control groups (Figure 8D). The correlations of the 22 immune cells with the six m6A-Ferr-related signature genes were calculated and individually plotted in lollipop plots (Figures 8E–J). According to p_adj <0.5 and |cor| >0.3, naive B cells, CD8+ T cells, regulatory T cells (Tregs), and activated natural killer (NK) cells were positively correlated with CDO1 and NOX4 but negatively correlated with ATG7, CYBB, and SLC3A2. In addition, AGPAT3 was negatively correlated with naive B cells, CD8+ T cells, and activated NK cells, and SLC3A2 was positively correlated with activated CD4+ memory T cells.