Identification of Potential Ferroptosis-Related Biomarkers and Immune Infiltration in Human Coronary Artery Atherosclerosis.

Objective: Ferroptosis is a specific subtype of programmed cell death, which plays an essential role in the immune-associated disease, atherosclerosis (AS). The purpose of this study was to identify potential ferroptosis-related gene biomarkers and its association with immune infiltration characteristics in atherosclerosis with bioinformatics methods.
Methods: Differentially expressed genes (DEGs) between AS and control groups were screened from GSE40231, analyzed for functional enrichment and then intersected with ferroptosis-related genes. Then, a random forest model was constructed based on these differentially expressed ferroptosis-related genes (DE-FRGs) and validated with dataset GSE132651. The performance of the models was evaluated with the area under receiver operating characteristic curves (AUC). Finally, we analyzed the correlation between DE-FRGs above and the characteristics of immune infiltration via CIBERSORT method.
Results: Six DE-FRGs (IL6, ANGPTL7, CDKN1A, AKR1C3, NOX4 and VLDLR) were detected based on dataset of GSE40231. Furthermore, a random forest model was constructed based on them with a compelling diagnostic performance of AUC = 0.8974 in the validation dataset GSE132651. In addition, the proportion of follicular helper T (Tfh) cells was significantly higher in AS group (P < 0.001). And we found significant correlation relationship between Tfh and expression level of ANGPTL7 (R = 0.35, P < 0.01), CDKN1A (R = 0.4, P < 0.0001), AKR1C3 (R = 0.64, P < 0.0001), NOX4 (R = 0.32, P < 0.01) and VLDLR (R = -0.43, P < 0.0001).
Conclusion: This study identified 6 DE-FRGs and validated a predicted model for the early prediction of AS, which also proved the close relationship between ferroptosis and immunity in the pathogenesis of AS.

Introduction

Coronary artery disease (CAD) resulting from AS is a leading cause of mortality worldwide. In 2017, it was estimated that CAD affects around 126 million individuals all over the world, accounting for 1.72% of the global population.1 AS is characterized by endothelial dysfunction or death in the early stage, among which oxidative stress is an important initiating factor.2 In the end, large amounts of lipids and several types of immune cells, including macrophages, T and B lymphocytes accumulated and narrowed down the arterial lumen by forming the atherosclerotic plaques in the arterial wall.3,4 Experimental studies have demonstrated that large amounts of macrophages are activated by different cytokines and mediators and become foam cells in the lumen lesion.5 T cells, despite its minority proportion, also have a very significant function in the immune regulation of AS.6 In addition, all subsets of T cells (including CD4+, CD8+, NK and follicular helper T cells) have been identified in human atherosclerotic plaque.6,7 Inflammatory cells released specific cytokines and enzymes, which finally promote plaque erosion and rupture, leading to partial or total occlusion of the affected artery.8

Studies over the past several decades have defined and characterized several cell death types such as apoptosis, pyroptosis, necroptosis, and ferroptosis, which have all been reported to accelerate the process of atherosclerotic lesion formation.9–13 The roles of ferroptosis in AS are attracting increasingly attention at present. Ferroptosis is characterized by overwhelming, iron-dependent accumulation of lipid hydroperoxides. Glutathione peroxidase (GPX4) is an essential enzyme that defends cells against oxidative damage. Specifically, it could reverse lipid peroxidation and ferroptosis by consuming glutathione. GPX4 inactivation and the depletion of intracellular glutathione blocked the metabolic pathway of lipid oxides, which induced iron homeostatic limits and occurrence of severe oxidative damage, thereby resulting in the progression of ferroptosis.14 Nowadays, more and more literatures support an important role of ferroptosis in the pathogenesis of AS. For instance, Yang et al reported that overexpression of prenyl diphosphate synthase subunit 2 (PDSS2) suppressed the ferroptosis of human coronary artery endothelial cells by promoting the activation of nuclear factor erythroid 2-related factor-2 (Nrf2) pathways.15 Moreover, trapping ferroptosis by effective and specific small molecular ferroptosis inhibitors such as deferoxamine (DFO) could alleviate AS through attenuating lipid peroxidation and endothelial dysfunction.13,16

The underlying possible molecular mechanisms of ferroptosis and its association with immunity in the field of AS are still in the veil. To identify key ferroptosis-related genes, we used integrated bioinformatic analysis methods to screen differentially expressed genes (DEGs) between AS and control artery tissue of CAD patients. Meanwhile, functional enrichment analysis was exhibited to uncover biological functions and enriched pathways for these genes. Besides, six DE-FRGs were obtained by intersecting these DEGs with the ferroptosis gene dataset from FerrDb database. Additionally, we established a random forest predictive model based on these DE-FRGs and externally verified it using another data set (GSE132651). Pearson correlation coefficient was also calculated to exhibit the association between immune cells and DE-FRGs. In this study, systematic bioinformatics analyses are applied to screen hub ferroptosis-related genes and its relationship with immunity in AS. Our results will help to shed light on new thoughts in the mechanism of CAD.

Materials and Methods

Microarray Data Download and Preprocessing

We downloaded the dataset of GSE40231 from GEO database, which was original from the Stockholm Atherosclerosis Gene Expression (STAGE) study. The primary aim of the study was to illustrate functionally associated genes in CAD. In this gene expression study, samples of atherosclerotic and unaffected arterial wall (distal part of mammary artery which was used as a graft) were isolated from 40 CAD patients undergoing coronary artery bypass grafting (CABG) surgery.17 After the microarray raw data were downloaded successfully, they were normalized by using a robust multichip average (RMA) background correction algorithm provided by Affy and limma package. Details of the data analysis workflow for the study design were described (Figure 1).

Figure 1

Flow chart for the whole study.

Functional Enrichment Analysis

GO term enrichment analysis, which contains subsections of biological processes, cellular component and molecular functions,18 was conducted through EnrichR web tool ().19 Critical signaling pathways enriched in DEGs were identified from the database of KEGG (Kyoto Encyclopedia of Genes and Genomes)20 on the same platform. P < 0.05 was considered statistically significant for both GO and KEGG analysis. Top 10 listed pathways with the smallest P value are displayed.

Recognition of Differentially Expressed Genes

The expression profiles of AS and unaffected arterial wall were compared to screen the DEGs. The cut-off criteria is as following: |log2(fold-change)|>1 and adjusted P-value <0.05. Volcano plot and heat map were drawn to display the expression of DEGs between AS and control groups with the “limma” package.21

Identification of DE-FRGs

We obtained a gene dataset which included 269 genes from the Ferroptosis Database ().22 Intersection of the ferroptosis genes and DEGs from GSE40231 was exhibited using the VennDiagram package in R. Principal component analysis (PCA) was applied for dimensionality reduction to distinguish the AS and control groups using the ggfortify package in R. Finally, the ggpaired package in R was utilized to exhibit the expression of these DE-FRGs in different groups using a paired sample t-test as the statistical analysis method. The P-value <0.05 was considered statistically significant.

Construction of Random Forest Model

The random forest is a machine learning technique of training and predicting samples with high accuracy via constructing a multitude of decision trees.23 The prediction model has been widely used for identifying potential predictors with the randomForest package24 in R. To achieve less biased or less optimistic performance of the training set, K-fold cross-validation with cv = 10 was applied. In the ongoing study, six DE-FRGs and disease state of control and AS were viewed as the classification feature and variables, respectively. The samples in GSE40231 (n = 80) were randomly assigned into the nonoverlapping training set and test set with the ratio of 2:1. Suitable values for key parameters were specified, such as mtry and ntree. The model built by the training set was verified by the test set further. The error rate was calculated to evaluate the combined classification through Out-of-bag (OOB) algorithm.25 The mean decrease accuracy (MDA) and the mean decrease Gini (MDG) were positively correlated with the importance of variables.26 Finally, GSE132651 was used to externally verify the model and demonstrate the prediction ability of this diagnostic model. ROC curves were visualized through pROC package.

Immune Infiltration Characteristic and Its Correlation with DE-FRGs

Immune infiltration characteristics of 22 types of immune cells were assessed by the CIBERSORT algorithm,27 including macrophages, T cells, natural killer (NK) cells, mast cells, B cells, dendritic cells (DC), monocytes, plasma cells, neutrophils, basophils and eosinophils. After obtaining the expression proportion matrix of immune cell subtypes with the guide of the CIBERSORT website, we use “ggplot2” package to draw boxplot to depict the distribution of immune cells. The Pearson correlation coefficient between each kind of immune cells and DE-FRGs was calculated, respectively, and the results with correlation coefficient and P value were exhibited using “ggpubr” package in R software.

Result

Identification and Functional Analysis of DEGs

The 80 samples which contained the AS and control groups were normalized via RMA algorithm. Three hundred and sixteen DEGs were identified after setting the cutoff values (|log2(fold-change)|>1 and adjusted P-value <0.05), including 145 upregulated genes and 171 downregulated genes. They were displayed in a volcano plot (Figure 2A). Results have also revealed the top five upregulated genes (CARTPT, SCG2, CYTL1, SLC10A4, HAPLN1) and top five downregulated genes (HOXC6, EMX2, HOXA5, RBPMS2, HOXA9) in the AS group. In addition, these DEGs were exhibited in a heat map (Figure 2B).

Figure 2

Identification of differentially expressed genes (DEGs) and functional enrichment analysis. (A) The DEGs identified were displayed in a volcano plot. Up-regulated genes are marked in red, and down-regulated genes are labeled in blue. (B) Heat map of the DEGs. (C) Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis among these DEGs.

Functional enrichment analysis was carried out with Enrichr web tool. We analyzed GO terms and KEGG pathway enriched in 316 DEGs. A collection of GO terms and pathways according to -log10 (P value) was displayed through bar plot (Figure 2C). For biological pathway analysis, the most enrichment terms were regulation of angiogenesis and inflammatory response. Molecular function subsection data indicated that receptor ligand, cytokine and chemokine activity involved in the DEGs significantly. Cellular component study displayed significant enrichment of these genes in collagen-containing extracellular matrix. KEGG pathway analysis mapped these genes into viral protein interaction with cytokine and cytokine receptor.

Collection and Analysis of DE-FRGs

We collected a dataset of 259 genes after deduplication from the Ferroptosis Database (FerrDb), which included 108 drivers, 69 suppressors and 111 markers of ferroptosis. The common genes between ferroptosis genes from FerrDb and DEGs of GSE40231 were visualized through a Venn diagram. Six DE-FRGs were identified, including Interleukin 6 (IL6), Angiopoietin Like 7 (ANGPTL7), Cyclin Dependent Kinase Inhibitor 1A (CDKN1A), aldo-keto Reductase Family 1 Member C3 (AKR1C3), NADPH Oxidase 4 (NOX4) and Very Low Density Lipoprotein Receptor (VLDLR) (Figure 3A). The relevant information of these screened ferroptosis-related genes is elaborated in Table 1. A protein-protein-interaction (PPI) network is shown in Figure 3B. Results from PCA revealed that genes mentioned above could successfully distinguish AS patients from those healthy people (Figure 3C). All these DE-FRGs were exhibited to be greatly changed in the AS group in GSE40231. P-value <0.0001 was considered statistically significant (Figure 3D).

Figure 3

Analysis of the differentially expressed ferroptosis-related genes (DE-FRGs). (A) Overlapped DE-FRGs in Venn diagram; (B) protein-protein-interaction (PPI) network of 6 DE-FRGs. (C) principal component analysis of these DE-FRGs; (D) gene expression level of 6 DE-FRGs was compared in AS and control groups with paired t-test, respectively. Significance markers: ****P < 0.0001.

Table 1

DE-FRGs of Coronary Artery Atherosclerosis

Gene	Full Name	Role in Ferroptosis	logFC	P-value	adj.P.Val
IL6	Interleukin 6	Marker	2.092407	3.22E-12	1.19E-10
ANGPTL7	Angiopoietin Like 7	Marker	1.257568	2.92E-12	1.09E-10
CDKN1A	Cyclin Dependent Kinase Inhibitor 1A	Suppressor	1.018134	3.16E-13	1.42E-11
AKR1C3	Aldo-keto Reductase Family 1 Member C3	Suppressor	1.169341	1.80E-15	1.19E-13
NOX4	NADPH Oxidase 4	Driver	1.162017	2.94E-16	2.15E-14
VLDLR	Very Low Density Lipoprotein Receptor	Marker	−1.15079	2.22E-17	1.97E-15

Abbreviations: logFC, log fold change; adj.P.Val, adjusted P value.

The random forest regression classification model was constructed based on the 6 DE-FRGs. When mtry = 1, the false-positive rate was the lowest (about 8.35%) (Figure 4A). The error rate became stable when the number of decision trees reached about 200 (Figure 4B). The importance of variables was weighed according to two parameters, mean decrease accuracy and mean decrease Gini. It was found that AKR1C3 was the most important variables for AS prediction (Figure 4C). After we established the random forest model with train data, the testing data sets were verified to prove its high prediction cost with AUC of 0.996 (Figure 4D). Additionally, GSE132651 was selected to externally verify the predictive model with a compelling performance of AUC = 0.8974.

Figure 4

Construction and evaluation of random forest model based on the 6 DE-FRGs. (A) The scatter plot exhibited the association between the false-positive rate and mtry index (from 1 to 6). (B) Trend of the related errors according to the number of decision trees. (C) The relationship between variables and two parameters, mean decrease accuracy (the left panel) and mean decrease Gini (the right panel), respectively. (D) Evaluation of the prediction efficiency of the prediction model using ROC analysis in GSE132651.

Relationships Between DE-FRGs and Immune Infiltration

CIBERSORT algorithm was used to evaluate the immune infiltration characteristics of AS. The relative proportion of immune cell subtypes was displayed in the boxplot. As shown in Figure 5A, follicular helper T cells (Tfh) showed higher proportion significantly in AS group, compared with control group (P < 0.001). Meanwhile, there was significant difference in naïve B cells, naïve CD4+ T cells, gamma-delta T cells (P < 0.01) and resting CD4+ memory T cells (P < 0.05) between AS and control. Further analysis about Pearson correlation coefficient demonstrated close relationship between DE-FRGs and immune infiltration (Figure 5B). We found significant association between Tfh and expression level of ANGPTL7 (R = 0.35, P < 0.01), CDKN1A (R = 0.4, P < 0.0001), AKR1C3 (R = 0.64, P < 0.0001), NOX4 (R = 0.32, P < 0.01) and VLDLR (R = −0.43, P < 0.0001). Besides, naïve B cells were shown to be related to ANGPTL7 (R = 0.25, P = 0.033) and CDKN1A (R = 0.38, P = 0.0011). And naïve CD4+ T cells was related to IL6 (R = −0.48, P = 0.0086) and ANGPTL7 (R = −0.42, P = 0.0025). The correlation relationship between other immune cell subtypes and 6 DE-FRGs did not show significant association (which was not shown). Overall, these 6 DE-FRGs were in varying degrees related to immune cells, especially Tfh.

Figure 5

Identification of immune infiltration and its correlation relationship with DE-FRGs. (A) Proportion of immune cell subtypes exhibited in boxplot diagram with Wilcoxon test. Significance markers: *P < 0.05; **P < 0.01; ***P < 0.001. (B) Pearson correlation coefficient analysis between immune cells and DE-FRGs. P values of P < 0.05 were labeled with red color.

Discussion

The term ferroptosis was newly proposed in 2012, which was defined as a specific cell death type induced by the small molecule erastin.14 It has been reported to be closely involved in the pathophysiological process of diseases characterized by sterile inflammation, including cancers, degenerative brain diseases and cigarette smoke-induced chronic obstructive pulmonary disease.28 However, definite evidence for the underlying molecular mechanisms of ferroptosis in the field of atherosclerosis is still limited. In the current study, we used integrated bioinformatics methods to identify key ferroptosis-related genes in the pathogenesis of atherosclerosis. The six ferroptosis-related genes (IL6, ANGPTL7, CDKN1A, AKR1C3, NOX4 and VLDLR) were speculated to play a vital role in CAD.

A random forest model with 10-fold cross-validation was constructed base on these hub genes. Cross-validation is a resampling procedure used to evaluate the machine learning model when sample sizes are limited. It usually results in a less biased or less optimistic generalization performance.29 The random forest model constructed in this study was verified further in GSE132651, which was published by Hebbel et al.30 In GSE132651, blood outgrowth endothelial cells (BOEC) from subjects with abnormal and normal coronary endothelial function were used for gene expression analysis. Our study addresses the feasibility of applying the model above to predict endothelial function with a satisfying diagnostic performance of AUC = 0.8974. Since coronary endothelial dysfunction is the earliest clinically detectable form of atherosclerosis,31 our study may also have practical clinical implications in diagnosis of atherosclerosis at an early stage.

Previous literatures reviewed the relationship of these screened ferroptosis-related hub genes with atherosclerosis. Interleukin-6 (IL-6) is a widely used cytokine that causes inflammation and oxidative stress. A study in 1994 mentioned that IL-6 was expressed locally in coronary atherosclerotic plaques, with expression level of 10- to 40-fold in the atherosclerotic arteries compared with the non-atherosclerotic artery.32 Not only IL-6 receptor antibody but also IL-6 inhibitor has been proved to be useful to alleviate atherosclerotic lesion development.33 Li et al reported that inhibition of ANGPTL7 significantly reversed tumor necrosis factor-α induced oxidative stress and endothelial cell adhesion,34 which identified Angptl7 as a potential target in atherosclerosis. CDKN1A, also named p21, encodes a potent cyclin-dependent kinase inhibitor. However, the role of p21 in atherosclerosis is controversial. The study carried out by Merched and Chan revealed the proatherogenic effect of p21, since loss of p21 protected against atherosclerosis in apoE-/- mice.35 Some studies considered p21 to be as an antiatherogenic target that can be used to treat or prevent atherosclerosis and postangioplasty restenosis.36,37 AKR1C3 is one of aldo/keto reductase superfamily. While AKR1C3 is expressed mainly in human steroid-hormone-targeting tissues including prostate and lung, it is also constitutively expressed in artery endothelial cells.38 Moreover, AKR1C3 was involved in the process of reactive oxygen species (ROS) accumulation in endothelial cells.38 NOX4, which belongs to the NADPH oxidase family, is constitutively active in smooth muscle cells (SMCs) and endothelial cells.39 Tong et al reported that Nox4 plays a pro-atherosclerosis role by stimulating SMC proliferation, migration and inflammation.39 VLDLR belongs to LDL receptor family, and its association with atherosclerosis is complicated. On the one hand, VLDLR is involved in peripheral triglyceride uptake, exerting a proatherogenic effect. On the other hand, VLDLR in SMC likely has a protective role against atherosclerosis. Besides, both pro- and anti-atherogenic roles of VLDLR in macrophages have been reported.40

Accumulated evidences point to a functional role for Tfh cells in AS. Since the inducible T cell costimulatory ligand (ICOSL) signal provided by DCs or B cells is critical for the differentiation and maintenance of Tfh cells, αICOSL antibody can be used to downregulate Tfh.41 Gaddis et al reported that αICOSL antibody-treated group showed a 10-fold reduction in the numbers of Tfh cells derived from Treg in the aorta root of AopE−/− mice, accompanied by a 30% reduction in lesion area of atherosclerosis.42 Besides, atherosclerosis development was alleviated in Tfh cell-deficient mice (Bcl6fl/flcreCd4),42 which proved that Tfh cells are indeed directly pro-atherogenic, instead of through other immune cell population impacted by ICOS/ICOSL interactions indirectly. Moreover, it was shown that oxidized low-density lipoprotein (ox-LDL) contributed to plaque formation via facilitating the generation of Tfh cells, which were original from nTreg cells, via reduction of IL2Rα and enhancement of IL6R expression.42 However, reports about the relationship between ferroptosis and Tfh cells are limited, let alone in the field of AS. A newest related study was just published, which illustrated selenium-GPX4-ferroptosis axis, instead of apoptosis or pyroptosis, played a key role in regulating Tfh homeostasis in immunized mice and young adults after influenza vaccination.43 In our study, we found significant correlation relationship between Tfh and expression level of ANGPTL7, CDKN1A, AKR1C3, NOX4 and VLDLR. The direct evidences about the association were not mentioned in papers published previously. Further experiment research is needed to classify whether these ferroptosis-related genes and proteins participate in the pathogenesis of AS by regulating local immune cells in artery wall, especially Tfh.

In this study, we have constructed and validated a ferroptosis-related genes associated model in AS. Our results revealed the powerful prediction ability of ferroptosis gene biomarkers and reflected their relationship with the immune infiltration. Nevertheless, the shortcomings of this study should also be pointed out. Firstly, our results were dependent on existing data by bioinformatic analysis, suggesting that scant analytical methods may limit the predictive capability of the constructed model and new proof from future studies may renew current results. Secondly, our analyses are based on public datasets, and further researches about the detailed molecular mechanism are necessary.

Conclusion

In the current study, we screened six DE-FRGs (IL6, ANGPTL7, CDKN1A, AKR1C3, NOX4 and VLDLR) in artery wall between the AS and control groups. Furthermore, a random forest model of AS was developed with these hub genes and verified in dataset GSE132651 (AUC = 0.8974). Immune cell infiltration and its correlation relationship with DE-FRGs were recognized. The DE-FRGs associated model provided a simple tool for the early prediction of AS, which improved our understanding of the association of ferroptosis and immunity in the pathogenesis of AS further.