Identification of Ferroptosis-Related Hub Genes and Their Association with Immune Infiltration in Chronic Obstructive Pulmonary Disease by Bioinformatics Analysis

Introduction

Chronic obstructive pulmonary disease (COPD) has an increasing rate of incidence in recent years and causes three million deaths annually, which brings about a heavy economic burden.¹ Currently, there are no effective target drugs applied to clinical practice so it is urgent to mine promising drug targets. Airway inflammation is an important feature and contributes to the pathogenesis and progression of COPD.² Ferroptosis refers to the programmed cell death induced by lipid peroxidation via iron-dependent pathway with unique morphological and biological features.³ Usually under environmental stresses or intra/inter-cellular signaling, many metabolic products such as reactive oxygen species (ROS) and phospholipid containing polyunsaturated fatty acid chain(s) (PUFA-PL) can trigger phospholipid peroxidation.⁴ Previous study proved that ferroptosis was involved in the pathogenesis of COPD.⁵ Stimulated by cigarette smoke, bronchial epithelium produced reactive oxygen species, which induced lipid peroxidation, membrane damage and even ferroptosis.^6–8 Glutathione peroxidase 4 (GPX4) – a vital antioxidant regulator – is also impaired during ferroptosis.⁹ Cigarette smoke extract altered ferroptosis-related genes expression in bronchoalveolar epithelial cells. Hypermethylation of the nuclear factor erythroid 2-related factor 2 (Nrf2) promoter could inhibit Nrf2/GPX4 axis, thus affecting ferroptosis in COPD.¹⁰ Otherwise, many studies suggest that various immune cells play vital roles in chronic airway diseases such as COPD. Innate immune cells, which were enhanced in small airways, modulated airway inflammation and remodeling.¹¹ Previous study reported that CD8+ T cells enhanced ferroptosis-specific lipid peroxidation in tumor cells.¹² The proliferation of B cells and antibody production was influenced by iron ion regulating the expression of Cyclin E1.¹³ Macrophages recognized oxidized phospholipids on the cell surface to clear ferroptosis cells via toll-like receptor 2 (TLR2).¹⁴ However, the interaction between ferroptosis and immune cells infiltration in COPD pathogenesis remains unclear.

The aim of this research is to identify ferroptosis-related hub genes and their association with immune cells infiltration in COPD lung tissues compared with normal ones. Additionally, we intend to construct interactive networks of hub genes with miRNAs, transcription factors and signal molecules and evaluate the diagnostic values of hub genes.

Materials and Methods

Data Acquisition

The mRNA expression microarray data of GSE38974,¹⁵ including 23 patients with COPD and 9 normal controls, were extracted from the Gene Expression Omnibus (GEO) datasets.¹⁶ The platform was GPL4133 Agilent-014850 Whole Human Genome Microarray 4x44K G4112F (Feature Number version). Lung tissues from 9 smokers with no evidence of obstructive lung disease and 23 smokers with COPD were examined for mRNA expression. All the clinical information including age, gender, sample source, smoking history, GOLD stage and FVC group was publicly accessible in GEO database (https://www.ncbi.nlm.nih.gov/geo/geo2r/?acc=GSE38974&platform=GPL4133).

The 259 ferroptosis-related genes (FRGs) were downloaded from the FerrDb database.¹⁷

Identification of Differentially Expressed FRGs

The GEO2R is a web tool internally stalled in GEO database specialized for analyzing differentially expressed genes between experimental group and control group. The GEO2R analysis between COPD group and control group was performed on the GEO datasets and the result of differential expression analysis was downloaded for further analysis. The cut‐off criteria for differential gene expression were the absolute value of log fold change (FC) >1 and P value <0.05. The gene list of differential expression analysis and FRGs were intersected to obtain the differentially expressed FRGs. The expressions of differentially expressed FRGs were plotted using the R package Complex Heatmap.¹⁸ The expression differences of differentially expressed FRGs between COPD group and normal group in GSE38974 were compared using Kruskal–Wallis test and Dunn’s test. The correlation analysis of differentially expressed FRGs was performed using Spearman’s correlation.

Gene Ontology (GO) Terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathways Analyses

The GO and KEGG enrichment analyses were conducted using the R package clusterProfiler.¹⁹ The screening criteria for significant terms were adjusted P values less than 0.05 and q values less than 0.2. Combined with logFC values, the enrichment analyses were performed by calculating the Z-scores using the R package GOplot.²⁰

STRING Database and Cytoscape Software

The STRING database is an online tool for analyzing protein–protein interaction. The PPI analysis was carried out using the STRING database.²¹ Cytoscape is a computer software that graphically displays, analyzes and edits the network, which contains multiple plugins. The results were processed and visualized in Cytoscape software (version 3.9.0). The key module was screened by the Molecular Complex Detection (MCODE) plugin. The differentially expressed FRGs were ranked by degrees and the top five genes were considered to be hub genes by cytoHubba plugin.

The Comparative Toxicomics Database

The Comparative Toxicomics Database (CTD) provides integrated information on complex interactions among chemical exposures, genes, proteins and diseases.²² In this study, we used it to estimate the inference scores of hub genes in several respiratory tract diseases.

miRNet, NetworkAnalyst and Encyclopedia of RNA Interactomes (ENCORI)

There are multiple online tools to analyze the interaction between non-coding RNAs and genes and predict target molecules. The miRNet is a useful online tool centering around miRNAs and their interacting molecules.²³ The NetworkAnalyst is an integrated and powerful database for gene expression analysis and construction of interacting networks.²⁴ In this article, they were utilized to explore and visualize the networks between hub genes and miRNAs, transcription factors and signal molecules. The ENCORI focuses on predicting RNA interaction.²⁵ Here, it was used to predict the upstream molecules lncRNAs targeting screened miRNAs. The screening criterion was set as strict stringency (the number of Ago CLIP-seq experiments is no less than five) and the top three lncRNAs were selected.

The Receiver Operating Characteristic (ROC) Curves of Hub Genes

The logistic regression model of hub genes was constructed using glm function in R software, and the ROC curves were plotted using R package pROC. The ROC curve of each hub gene can help us determine whether its expression has diagnostic value to some extent. The true positive rate (TPR) or sensitivity refers to the number of true positive samples detected divided by the number of all true positive samples. The false-positive rate (FPR) refers to the number of false-positive samples detected divided by the number of all true negative samples. Specificity refers to the number of true negative samples detected divided by the number of all true negative samples. The abscissa represents 1 – specificity and the ordinate represents sensitivity. The area under the curve (AUC) is used to determine the prediction accuracy. The AUC is usually between 0.5 and 1.0. The ROC curve has low/moderate/high accuracy when the AUC is 0.5~0.7/0.7~0.9/more than 0.9, respectively. The Youden’s index (sensitivity plus specificity minus one) is used to assess the authenticity of the model. As it gets closer to 1.0, the model is much more authentic.

CIBERSORT Algorithm

CIBERSORT is an R/web tool for deconvolution of expression matrices of human immune cell subtypes based on the principle of linear support vector regression.²⁶ By the way of the CIBERSORT algorithm, we analyzed the proportions of 22 types of immune cells infiltration in patients with COPD and normal controls. The infiltration differences in patients with COPD and normal controls were compared using Kruskal–Wallis test and Dunn’s test. The Spearman correlation analysis was carried out to show the correlation within differentially infiltrated immune cells and the association between hub genes and differentially infiltrated immune cells.

Statistical Methods

All the statistical calculations were conducted in R software (version 3.6.3). The corresponding R packages were described as above. The statistical significance was marked as follows: ns, p≥0.05; *p< 0.05; **p < 0.01; ***p < 0.001. A p.adjust was the corrected p value obtained by the p value correction method; a q-value was an adjusted p-value, taking into account the false discovery rate (FDR). As the p value/p.adjust/q-value is less than 0.0001, scientific notation (exponent, E) is used. For instance, 0.0000267 is written as 2.67E-05.

Results

Identification of Differentially Expressed FRGs

The design of this research was shown in the flow chart (Figure 1). Principal component analysis (PCA) was conducted to show that there was a good degree of clustering between the two groups (Figure 2A). After intersection, 102 genes were obtained (Figure 2B). Under the condition of absolute values of logFC > 1 and P values <0.05, 15 differentially expressed FRGs were discovered including 11 upregulated genes and 4 downregulated genes (Table 1). The 15 differentially expressed FRGs between COPD and normal groups were presented in volcano plot and heatmap (Figure 2C and D). The volcano plot showed the distribution of gene expression between COPD and normal groups. Genes with an adjusted P-value <0.05 and absolute fold-change value > 1 were considered as differentially expressed genes. Each point represented one gene. Red dots indicated significantly upregulated genes and blue dots indicated significantly downregulated genes. The top three upregulated genes included IL6, ATM and TNFAIP3 and the top two downregulated genes included IL33 and TGFBR1. Furthermore, the expression differences of 15 differentially expressed FRGs between the COPD and normal groups in GSE38974 were shown in box plots (Figure 3A). As shown in Figure 3A, 15 ferroptosis-related genes were all significantly differentially expressed between COPD group and normal samples, which was consistent with Figure 2C and D. To explore the expression correlation of these ferroptosis-related genes, correlation analysis was performed. The Spearman correlation analyses of 15 differentially expressed FRGs were shown in heatmap (Figure 3B). As shown in Figure 3B, the expressions of most up-regulating genes were significantly positively correlated with each other, so were the down-regulating genes. Similarly, the expressions of the most up-regulated genes had a significantly negative correlation with those of the four down-regulated genes.

Table 1 The Differentially Expressed FRGs in COPD Group Compared with Normal Group

Figure 1 Flow chart.

Figure 2 Identification of differentially expressed ferroptosis-related genes (FRGs). (A) The principal component analysis (PCA) plot of samples in GSE38974. (B) Venn diagram of GEO2R result and FRGs. (C) Volcano plot of differential genes between COPD and normal groups. The top five genes were labelled (upregulated-IL6, ATM and TNFAIP3, downregulated-IL33 and TNFBR1). (D) Heatmap of 15 differentially expressed FRGs.

Figure 3 Expression and correlation of differentially expressed FRGs. (A) Box plot of 15 differentially expressed FRGs in COPD group compared with normal group. The significance markers are shown as: *, P<0.05; **, P<0.01; ***, P<0.001. (B) Heatmap of correlation of 15 differentially expressed FRGs.

The GO and KEGG Enrichment Analyses of Differentially Expressed FRGs

To analyze the potential biological functions of these differentially expressed ferroptosis-related genes, we carried out GO and KEGG enrichment analyses by way of R software. In total, 739 biological processes (BPs), 11 cellular components (CCs), 21 molecular functions (MFs) and 26 KEGG pathways were enriched. The GO term results exhibited that differentially expressed FRGs were mainly involved in regulation of smooth muscle cell proliferation, membrane microdomain, membrane raft, caveola, cytokine receptor binding, cytokine activity, and transforming growth factor beta receptor binding. The KEGG analysis indicated that differentially expressed FRGs participated in cell senescence pathway, FoxO signaling pathway and HIF1 signaling pathway (Figure 4A and B, Table 2). After combining with expression levels (logFC), the Z-scores showed that all the significant terms could be positively regulated by differentially expressed FRGs (Z-scores >0) (Figure 4C and D). These findings implied that 15 differentially expressed FRGs may participate in inflammatory responses and airway remodeling in COPD pathogenesis.

Table 2 The Most Significant Terms of GO and KEGG Enrichment Analyses

Figure 4 GO and KEGG enrichment analyses of differentially expressed FRGs. (A) Lollipop plot of significant terms. (B) Circular network of significant terms and genes. Blue nodes represent terms, red nodes represent genes, and connecting lines represent the relationship between terms and genes. (C) Bubble plot of significant terms combined with logFC. A Z-score greater than zero indicates positive regulation, a Z-score less than zero indicates negative regulation, and absolute value of Z-score represents the probability of regulation. (D) Donut plot of significant terms combined with logFC. Each column of the inner circle corresponds to one term, the height of column represents adjusted P value, and the filled color represents Z-score of each term.

Construction of Protein–Protein Interaction (PPI) Network and Identification of Key Module and Hub Genes

To determine the interactive relationship among differentially expressed FRGs, the protein–protein interaction analysis was conducted. The interaction of 15 candidate genes was analyzed in STRING database, and the results were visualized in Cytoscape software. The results showed that these differentially expressed FRGs interacted with each other (Figure 5A) and displayed the interaction number of each gene (Figure 5B, Supplementary Table 1). In total, there existed 15 nodes and 134 edges. The MCODE plugin analysis showed that there existed one key module containing 11 nodes and 50 edges including GDF15, IL6, ATF3, PTGS2, TGFBR1, HIF1A, CDKN1A, ATM, HMOX1, TNFAIP3 and MYB (Figure 5C). The cytoHubba plugin analysis identified five hub genes, including HIF1A, IL6, PTGS2, CDKN1A and ATM (Figure 5D). The Venn diagram indicated the overlap of predicted hub genes (Figure 5E). The detailed information of five hub genes can be seen in Supplementary Table 2.

Figure 5 PPI network, key module and hub genes of differentially expressed FRGs. (A) The PPI among 15 differentially expressed FRGs. (B) The interaction number of each differentially expressed FRG. (C) Key module of the PPI network screened by MCODE plugin. (D) Hub genes screened by cytoHubba plugin. (E) The overlap of predicted hub genes.

Evaluation of Correlation Between Hub Genes and Respiratory Tract Diseases

In order to estimate the theoretical association between predicted hub genes and chemical/environmental exposures, the five hub genes were analyzed in the Comparative Toxicomics Database and four respiratory tract diseases were chosen including COPD, chronic bronchitis, pulmonary emphysema and non-small cell lung cancer (NSCLC). The average inference scores of five hub genes in COPD (46.85) were higher than those in chronic bronchitis (35.72) and pulmonary emphysema (17.17) but lower than those in NSCLC (55.56) (Figure 6). The findings implied that five predicted hub genes might participate in multiple pathophysiological processes in respiratory diseases.

Figure 6 The correlations between hub genes and respiratory tract diseases in comparative toxicomics database.

Construction of the Networks Between Hub Genes with miRNAs, Transcription Factors and Signal Molecules

The hub genes could probably play a role by acting as transcription factors and vital signal molecules or interacting with intracellular non-coding RNAs. In order to predict upstream or downstream molecules of five hub genes and speculate on the mechanism of action of each hub gene, interactive network analysis was conducted. The five hub genes were uploaded to the miRNet online database to analyze the interaction with miRNAs and transcription factors in human lung tissues. In total, 44 miRNAs were predicted and two miRNAs, hsa-let-7b-5p and hsa-miR-1-3p, both targeting five hub genes, were selected for further exploration (Figure 7A, Supplementary Table 3). The transcription factor-hub gene regulatory network consisted of 217 interactions between 164 transcription factors and five hub genes. Five transcription factors including EGR1, NFKB1, RELA, SP1 and STAT3, which had the highest connectivity with hub genes, were selected (Figure 7B, Supplementary Table 3). The ENCORI database was used to screen upstream lncRNAs of the two miRNAs. Six lncRNAs were predicted: NUTM2A-AS19, XIST and NEAT1 (targeting hsa-let-7b-5p); RMRP, MALAT1 and AL162431.2 (targeting hsa-miR-1-3p) (Figure 7C, Supplementary Table 3). Next, five hub genes were uploaded into the NetworkAnalyst database to analyze the interaction with signal molecules. TP53 was prominent due to interacting with three hub genes in a network of signal molecules (Figure 7D, Supplementary Table 3).

Figure 7 The interaction network of hub genes in miRNet and network Analyst. (A) The network of hub genes with miRNAs. (B) The network of hub genes with transcription factors. The fuchsia nodes represent hub genes and the green nodes represent transcription factors. The five transcription factors that connect with at least four hub genes are labelled. (C) The network of hub genes with signal molecules. The signal molecules that connect with at least two hub genes are labelled. (D) The predicted lncRNA-miRNA-hub gene regulatory network. Yellow diamonds represent lncRNAs, green ellipses represent miRNAs, and blue rectangles represent hub genes.

The ROC Curves of Hub Genes

To determine the diagnostic value in discriminating COPD patients from normal controls, the ROC curves of each hub gene were plotted using R software. The logistic regression model of hub genes was constructed based on glm function. The formula was “-88.166 + 3.0089*HIF1A + −2.8988*IL6 + 2.8957*PTGS2 + 3.2435*CDKN1A + 7.3934*ATM”. As shown in Figure 8, the expression levels of HIF1A (AUC: 0.923, CI: 0.804-1.00) and ATM (AUC: 0.976, CI: 0.926-1.000) had high predictive accuracy (Figure 8A and E). The expression levels of IL6 (AUC: 0.826, CI: 0.608-1.000) and CDKN1A (AUC: 0.860, CI: 0.653-1.000) had moderate predictive accuracy (Figure 8B and D). The expression level of PTGS2 had low predictive accuracy (AUC: 0.681, CI: 0.471-0.892) (Figure 8C). The AUC of combination of five hub genes was 0.981 (CI: 0.940-1.000) (Figure 8F). When the cut-off threshold was 1.398, the sensitivity, specificity and Youden index were 0.957, 1.000 and 0.957, respectively. These results indicated that this model had high accuracy and authenticity to distinguish COPD group from normal group.

Figure 8 The receiver operating characteristic (ROC) curves of hub genes. (A) ROC curve of HIF1A. (B) ROC curve of IL6. (C) ROC curve of PTGS2. (D) ROC curve of CDKN1A. (E) ROC curve of ATM. (F) ROC curve of five genes combination.

The Immune Cells Infiltration Characteristics in Patients with COPD and Normal Controls

The infiltrating status of various immune cells in lung tissues had obvious differences (Figure 9A). Monocytes and macrophages accounted for the majority of all infiltrating cells, especially in COPD lung tissues. The infiltration differences in both groups are shown in Figure 9B. Seven types of immune cells, including CD8 T cells, activated NK cells, monocytes, M0 macrophages, M2 macrophages, resting dendritic cells and resting mast cells, had differential infiltration in patients with COPD compared with normal controls. The adjusted P-values of seven kinds of immune cells were 0.002, 0.001, 0.025, <0.001, 0.002, 0.008 and 0.020, respectively. Monocytes and M0 macrophages were upregulated in COPD lung tissues, while CD8 T cells, activated NK cells, M2 macrophages, resting dendritic cells and resting mast cells were downregulated. Figure 10A reveals the correlations between differentially infiltrated immune cells. Monocytes had positive correlations with CD8 T cells, M2 macrophages and resting dendritic cells (r=−0.39, −0.67 and −0.46, respectively). M0 macrophages had inverse correlations with CD8 T cells and activated NK cells (r=−0.38 and −0.58, respectively). However, CD8 T cells had positive correlations with M2 macrophages, resting dendritic cells and resting mast cells (r = 0.54, 0.40 and 0.56, respectively). Resting mast cells were positively associated with M2 macrophages and resting dendritic cells (r = 0.36 and 0.61, respectively). The correlations between the expression of hub genes and differentially infiltrated immune cells are displayed in Figure 10B. Positive associations were observed between monocytes and IL6, monocytes and PTGS2, monocytes and CDKN1A (r = 0.53, 0.42 and 0.62, respectively). M0 macrophages were also positively associated with HIF1A and ATM (r = 0.50 and 0.52, respectively). However, CD8 T cells were strongly negatively associated with HIF1A, IL6, PTGS2 and CDKN1A (r=−0.68, −0.83, −0.72 and −0.79, respectively). The remaining several types of immune cells also had weakly to moderately negative correlations with the expression of most of the hub genes as displayed in the heatmap.

Figure 9 Immune infiltration of COPD lung tissues compared with normal tissues. (A) Stack bar chart of proportions of the immune cells infiltration. (B) Box plot of proportions of the immune cells infiltration. The significance markers are shown as: ns, P>0.05; *, P<0.05; **, P<0.01; ***, P<0.001.

Figure 10 Differentially infiltrated immune cells and hub genes. (A) Heatmap of correlations of differentially infiltrated immune cells. (B) Heatmap of correlations of hub genes with differentially infiltrated immune cells. The significance markers are shown as: *, P<0.05; **, P<0.01.

Discussion

Accumulating evidence indicates that ferroptosis participates in the pathogenesis of COPD. Previous review summarized that ferroptosis can affect inflammation through immunogenicity and ferroptosis inhibitors may benefit certain diseases through their anti-inflammatory effects.² However, more research is required to better our understanding of ferroptosis in pathogenesis of COPD. In our study, we obtained 15 differentially expressed FRGs in patients with COPD compared with normal controls through bioinformatics analysis. Several hub genes were reported in the previous study. For instance, HIF1A, as a switch gene, was upregulated in COPD cases using network-based analysis implemented by SWIM software.²⁷ CDKN1A played important functions in the development and progression of COPD.²⁸ In our study, the enrichment analyses of 15 differentially expressed FRGs were conducted to explore their potential functions. The results indicated that they were associated with airway inflammatory response and remodeling. For example, cell senescence occurs in many pathological processes in COPD, which is consistent with previous reports. Cell senescence impedes iron-mediated cell death pathways by impairing ferritinophagy, a lysosomal process that promotes ferritin degradation.²⁹

Next, we constructed a PPI network of 15 differentially expressed FRGs and first identified five ferroptosis-related hub genes, including HIF1A, IL6, PTGS2, CDKN1A and ATM. To further explore the correlation between hub genes and diseases, we analyzed the inference scores for four respiratory tract diseases in CTD and found that five hub genes were closely correlated with COPD and other respiratory tract diseases. These findings reminded us that it was vital to clarify the mechanism of action of these genes in COPD pathogenesis.

Bioinformatics methods provide us with a convenient way to predict crosstalk networks and screen potential biomarkers in COPD. A large number of miRNAs and lncRNAs were reported to be involved in COPD initiation and development. MALAT1/miR-146a/COX2 (namely PTGS2) axis affected the lung function of patients with COPD.³⁰ Some non-coding RNA targets including miR-195, miR-181c and TUG1 are viable for alleviating COPD in vivo.³¹ However, to our knowledge, previous articles reporting the correlation between non-coding RNAs and ferroptosis mainly focused on multiple cancers. The present study constructed the networks of hub genes with miRNAs and transcription factors in miRNet database and identified two key miRNAs, namely, hsa-let-7b-5p and hsa-miR-1-3p. The hsa-let-7b-5p participated in endothelial mitochondrial dynamics and acted as a biomarker for diagnosing Parkinsonian Syndromes.^32,33 The hsa-miR-1-3p inhibited lung adenocarcinoma cell tumorigenesis and improved gefitinib resistance in EGFR mutant lung cancer cell.^34,35 Additionally, the upstream lncRNAs of two miRNAs were predicted using ENCORI database and we found six lncRNAs with the most experimental evidence. Among them, NUTM2A-AS19 and AL162431.2 are newly reported in this study. XIST and MALAT1 played an important role in mitochondrial dysfunction, cell senescence and epigenetic alterations in COPD pathogenesis under the condition of tobacco smoke exposure.³⁶ NEAT1 promoted activation of inflammasomes in macrophages.³⁷ RMRP promoted the progression of NSCLC via competing with miR-1-3p.³⁸ Thus, we speculate that the following axes may regulate ferroptosis in COPD pathogenesis including NUTM2A-AS19 or XIST or NEAT1/hsa-let-7b-5p/hub gene axes and RMRP or MALAT1 or AL162431.2/hsa-miR-1-3p/hub gene axes. Moreover, the identification of five important transcription factors including EGR1, NFKB1, RELA, SP1 and STAT3 would be the groundwork for molecular mechanisms of ferroptosis in COPD pathogenesis. EGR1 was indispensable for MUC5AC expression induced by cigarette smoke in human bronchial epithelial cells.³⁹ Genetic knockdown of RELA (NFKB subunit) diminished IL6 production in HBE cells.⁴⁰ SP1 was crucial for anti-inflammatory molecule IL10 secretion in the phototherapy effect in HBE cells.⁴¹ STAT3 was a vital molecule in regulating the expression of inflammatory cytokines in COPD murine model.⁴² Notably, the signal molecule network revealed that TP53 connected with IL6, CDKN1A and ATM, suggesting that TP53 may be a potential driver of COPD towards lung cancer. This indicated that ferroptosis may also participate in COPD-related carcinogenesis.

Previous studies concentrated on the construction of a ferroptosis-related gene model for prognosis in cancer. For example, researchers screened ten ferroptosis-related genes, which served as potential prognostic biomarkers.⁴³ To testify the diagnostic values of hub genes, we conducted ROC analyses and discovered that each of them varied in predictive accuracy, while combination of five genes could serve as a fine model to distinguish patients with COPD from normal controls (AUC: 0.981, CI: 0.940-1.000).

Although immune infiltration in malignancies keeps attracting the attention of researchers, very few reports explored the immune infiltration in COPD. The spatially confined eosinophil-rich type 2 microenvironments were identified in COPD.⁴⁴ The proportion of T cells decreased in the lungs of current smokers and patients with COPD, whereas the proportion of macrophages increased.⁴⁵ In our study, we uncovered the immune infiltration status in patients with COPD compared with normal controls. Monocytes were the majority of immune cells in both groups and increased prominently in COPD group. Monocytes, as an essential part of innate immune system, influence human diseases both by direct effects and by differentiating into macrophages.⁴⁶ The cytokine response of monocytes to bacteria was compromised in smoking-induced COPD and thus impaired immune response.⁴⁷ Macrophages are plastic in response to various tissue microenvironment and external stimuli. We found that the proportion of M0 macrophages increased markedly, which could serve as reserves ready for polarization stimuli. M1 macrophages primarily take part in pro-inflammatory responses, however, they were not observed to increase remarkably in this study. M2 macrophages, which primarily participate in anti-inflammatory responses, decreased dramatically in COPD group, suggesting that their functions may be undermined in COPD pathogenesis. CD8 T cells were observed to decrease drastically, which was inconsistent with previous researches. The number of IFN-γ-producing CD8+ and CD4+ lymphocytes increased in the lungs of patients with COPD.⁴⁸ The frequencies of CD8+ T cell subsets increased observably in patients with COPD compared with normal controls and non-smokers.⁴⁹ It may be partial due to the difference between statistics-based bioinformatics methods and flow cytometry assays. Activated NK cells decreased in COPD lung tissues. However, evidence from other studies revealed that the proportions of NK cells increased in BAL fluid of patients with COPD.⁵⁰ Another two researches claimed that the number of NK cells in the lung parenchyma of patients with COPD was at the same level as that in the peripheral blood, and bronchoalveolar lavage fluid in healthy smokers.^51,52 It seems that NK cells in lung tissues may have different effects compared with those from blood and BALF. Moreover, we found resting dendritic cells and resting mast cells, which were the minority of infiltrating immune cells, decreased in patients with COPD. Dendritic cells present antigens and activate naive T and B cells.^53,54 Mast cells can interact with multiple immune cells and structural cells and thereby facilitate inflammatory responses, airway remodeling and angiogenesis.^55,56 The functions of them might be antagonized by monocytes and M0 macrophages in some degree in COPD lung tissues. The proportions of monocytes and M0 macrophages were positively associated with most of hub genes, whereas CD8 T cells, activated NK cells, M2 macrophages, resting dendritic cells and resting mast cells were negatively associated most of the hub genes. It can be inferred that the functions of monocyte and M0 macrophages may be promoted by these hub genes. The other types of infiltrating cells, for instance, CD8 T cells, were likely to be inhibited by these hub genes. HIF1A was validated to drive ferroptosis in clear cell carcinomas and ATM was essential for promoting ferroptosis.^57,58 IL6 and PTGS2 were confirmed as the downstream markers of ferroptosis.^59,60 CDKN1A was required to suppress ferroptosis.⁶¹ We speculate that these high-expressed hub genes in COPD group may get involved in the regulation of ferroptosis in structural cells of pulmonary parenchyma and thus affect the infiltrating immune cells residing in pulmonary interstitium or recruited from peripheral blood, which could lead to differentially histopathological changes in the lungs of patients with COPD. Taken together, the immune cells infiltration contributes to the pathogenesis of COPD in a sophisticated manner and more research is in urgent need to elucidate the situation.

Our study had obvious limitations. The number of cases included in our study was relatively small. Due to the lack of detailed clinical information, correlations between hub genes and clinical characteristics cannot be explored. Another apparent deficiency is that we did not perform basic experiments to validate the expression of hub genes and their correlation with immune cells. For now, our study can provide a theoretical basis for further explorations of ferroptosis-related phenotypes in COPD research.

Conclusion

We identified five ferroptosis-related hub genes (HIF1A, IL6, PTGS2, CDKN1A and ATM) in COPD, a combination of which had diagnostic value. Two miRNAs, five transcription factors and one signal molecule were predicted to target these hub genes, and the lncRNA-miRNA-hub gene regulatory network was constructed. Ferroptosis-related hub genes were significantly associated with immune infiltration in the lung tissues of patients with COPD.