{'fraction_non_alphanumeric': 0.028458349038823973, 'fraction_numerical': 0.021522804372408593, 'mean_word_length': 5.1872667910447765, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Identification of Ferroptosis-Related Hub Genes and Their Association with Immune Infiltration in Chronic Obstructive Pulmonary Disease by Bioinformatics Analysis

Purpose Ferroptosis and immune infiltration are involved in the pathogenesis of chronic obstructive pulmonary disease (COPD). We aim to identify ferroptosis-related hub genes and analyze their association with immune infiltration in COPD through bioinformatics methods. Materials and Methods The mRNA microarray data of GSE38974 were downloaded from Gene Expression Omnibus to obtain differentially expressed genes (DEGs). The DEGs were intersected with ferroptosis-related genes (FRGs) from FerrDb to obtain differentially expressed FRGs. GO and KEGG enrichment and protein–protein interaction (PPI) analyses of differentially expressed FRGs were conducted in R software and STRING database. The key module and hub genes were screened by Cytoscape software. MiRNAs, transcription factors and signal molecules were predicted in miRNet and NetworkAnalyst. The disease correlation in the Comparative Toxicomics Database (CTD) and the receiver operating characteristic (ROC) curves of hub genes were analyzed. Immune infiltration was evaluated by CIBERSORT algorithm. Spearman correlation analyses were conducted between hub genes and differentially infiltrated immune cells. Results Fifteen differentially expressed FRGs were identified, which were enriched in some terms involving airway inflammatory responses and structural remodeling. Five hub genes were selected including HIF1A, IL6, PTGS2, CDKN1A and ATM. Inference scores in CTD indicated their association with COPD. Two miRNAs, five transcription factors and one signal molecule were predicted. The combination of hub genes could be a fine diagnostic indicator of COPD (AUC: 0.981, CI: 0.940-1.000). Immune infiltration evaluation showed that 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. The hub genes were significantly associated with differentially infiltrated immune cells. Conclusion We identified five ferroptosis-related hub genes (HIF1A, IL6, PTGS2, CDKN1A and ATM) in COPD, and found that they may influence the pathogenesis of COPD by regulating ferroptosis and thus affecting infiltrating immune cells.


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. 1 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. 2 Ferroptosis refers to the programmed cell death induced

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. 18 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.

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. 21 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. 22 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. 23 The NetworkAnalyst is an integrated and powerful database for gene expression analysis and construction of interacting networks. 24 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. 25 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 falsepositive 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. 26 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.

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  positively correlated with each other, so were the down-regulating genes. Similarly, the expressions of the most upregulated genes had a significantly negative correlation with those of the four down-regulated genes.

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.

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.

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.

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   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).

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  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.

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 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

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. 2 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. 27 CDKN1A played important functions in the development and progression of COPD. 28 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. 29 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. 30 Some non-coding RNA targets including miR-195, miR-181c and TUG1 are viable for alleviating COPD in vivo. 31 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. 36 NEAT1 promoted activation of inflammasomes in macrophages. 37 RMRP promoted the progression of NSCLC via competing with miR-1-3p. 38 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. 39 Genetic knockdown of RELA (NFKB subunit) diminished IL6 production in HBE cells. 40 SP1 was crucial for anti-inflammatory molecule IL10 secretion in the phototherapy effect in HBE cells. 41 STAT3 was a vital molecule in regulating the expression of inflammatory cytokines in COPD murine model. 42 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. 43 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. 44 The proportion of T cells decreased in the lungs of current smokers and patients with COPD, whereas the proportion of macrophages increased. 45 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. 46 The cytokine response of monocytes to bacteria was compromised in smoking-induced COPD and thus impaired immune response. 47 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. 48 The frequencies of CD8+ T cell subsets increased observably in patients with COPD compared with normal controls and non-smokers. 49 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. 50 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. 61 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.

1233
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 ferroptosisrelated 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.

Disclosure
The authors report no conflicts of interest in this work.