Identification of the Key Immune-Related Genes in Chronic Obstructive Pulmonary Disease Based on Immune Infiltration Analysis

Introduction

Chronic obstructive pulmonary disease (COPD) is among the top three diseases in the world with respect to morbidity and mortality. It has become an important public health problem because of its high morbidity, high mortality and heavy socioeconomic burden.¹ Although great progress has been made in the treatment of COPD, the existing treatment methods cannot block disease progression or reduce mortality. Thus, we need to explore new diagnostic biomarkers for COPD and its potential molecular mechanisms.

With the development of bioinformatics, recent studies have shown a correlation between susceptibility to COPD and gene expression.² For example, MicroRNA-23a may be a potential biomarker for differentiating frequent exacerbators from infrequent exacerbators of COPD.³ In addition, IL6 and SOCS3 are also considered to play a key role in COPD and could be used as therapeutic targets for COPD.⁴ Therefore, bioinformatic analysis of gene expression profiles may help to screen for new COPD biomarkers and thus improve the treatment of COPD.

Smoking is the main cause of COPD, and exposure to environmental particulate matter (PM) will also greatly increase the burden of COPD.⁵ Research has reported that cigarette smoke induces oxidative stress,⁶ which damages biomacromolecules, leading to cell dysfunction or cell death. Moreover, oxidative stress also destroys the extracellular matrix, causing an imbalance of proteases and antiproteases and promoting the inflammatory response, ultimately leading to COPD. The airway and lung tissue of patients with COPD triggers continuous congenital and adaptive immune inflammatory responses in which neutrophils, macrophages, T lymphocytes and other inflammatory cells are involved in the pathogenesis.^7,8 Infection, inflammation, cigarette smoke and immune injury products can induce lung autoimmunity, accelerating the pathogenesis of COPD. However, the specific immune pathway is still unclear. Traditional research methods for assessing immune cells in the past have included immunohistochemistry and flow cytometry, both of which inevitably have limitations when it comes to fully analysing the composition of immune cells. In addition, flow cytometry may lead to cytolysis of certain cell types due to the high sample handling requirements.^9,10 However, the immune response in tumours involves a large number of specific cell types.¹¹ To better understand the diversity and nature of infiltrating immune cells in COPD, enumerating the number of immune cells in an aggregated fashion is a prerequisite. Cibersort is a general purpose gene expression based deconvolution algorithm that quantifies cellular components from the gene expression profile of a tissue.¹² As a result, different types of infiltrating immune cells can be quantified simultaneously, allowing the method to avoid the concerns of various surface markers and cell separation. In this study, we identified differential immune cell infiltration between normal and COPD samples by using CIBERSORT. Based on this information, we determined the differential immune-related genes associated with differential immune cells as hub genes of COPD by correlation analysis. In addition, a logistic regression model was used to evaluate the diagnostic value of hub genes. Finally, we further explored the regulatory mechanism of transcription factors on hub genes. Therefore, the integrated analysis of immune cells and immune-related genes will help to identify diagnostic markers of COPD, thus contributing to the clinical treatment of COPD.

Materials and Methods

Data Collection

Matrix files of the GSE57148 dataset (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE57148) based on the GPL11154 platform were extracted from the Gene Expression Omnibus (GEO) database, which contained 98 COPD and 91 normal individuals.¹³ Besides, we downloaded the GSE76925 dataset (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE76925), which contains RNA sequencing data from 111 COPD cases and 40 normal lung tissues,¹⁴ for diagnostic marker expression assessment and validation of diagnostic value. Additionally, 1509 IRGs were obtained from the ImmPort database. This study was certified by the Medical Ethics Committee of Yan’an Hospital, Kunming, Yunnan Province, China: Since this study uses human COPD-related data from a public database (GEO database) and does not involve any human medicine studies, this study is not subject to ethical review.

Evaluation of Immune Cell Subtype Distribution

CIBERSORT is an analytic tool that transforms the normalized gene expression matrix into the composition of 22 immune cell types based on the deconvolution algorithm.⁹ In the present study, the composition of immune cells between COPD and normal samples was calculated with CIBERSORT. The algorithm was performed with the LM22 signature and 1000 permutations. Samples of the present study were accurately satisfied with P < 0.05, and 90 COPD samples and 90 normal samples were selected for our further analysis. The fractions of the 22 immune cell types evaluated added up to 1 for each sample.

Identification of Differentially Expressed Immune-Related Genes (DE-IRGs) in COPD

The gene expression matrix of the GSE57148 dataset was analyzed to obtain DEGs between the COPD patient samples and the normal samples with the “limma” package of R. Then, significance analysis of microarrays was further conducted, and |log₂ fold change (FC)| > 0.5 and false discovery rate (FDR) < 0.05 were set as the selection criteria for DEGs. We combined the DEGs with 1509 IRGs to obtain the overlapping genes, which were used for subsequent analysis.

Functional Enrichment Analysis and Correlation Analysis

The DE-IRGs were annotated according to the Gene Ontology (GO)¹⁵ function to explore the biological significance of each gene. The GO annotation mainly contained three parts: biological process (BP), molecular function (MF), and cellular component (CC). We then investigated the signaling pathways of these DE-IRGs on the basis of the Kyoto Encyclopedia of Genes and Genomes (KEGG).¹⁶ GO term and KEGG pathway analyses were conducted using the clusterProfiler function in the R package. The cutoff criteria of P < 0.05 and FDR < 0.05 were defined as significant. Subsequently, we analyzed differential immune cells and DE-IRGs with the cutoff criteria of Pearson correlation coefficient (r) > 0.6 and P < 0.05. The DE-IRGs that had correlations with differential immune cells were selected as hub genes. Transcription factors of hub genes were predicted by human TF information (Network Analyst, http://www.networkanalyst.ca), and the regulatory network of transcription factors (TFs) and hub genes was visualized by Cytoscape software.

Ingenuity Pathway Analysis (IPA)

IPA is an integrated software based on cloud computing that can analyze, integrate and mine experimental data from gene expression, microRNA, SNP microarray, metabolomics and proteomics. The ingenuity pathway knowledge base (IPKB), as a key part of IPA, collects millions of pieces of information about the interactions between proteins, genes, tissues and diseases. We uploaded the list of DE-IRGs to IPA software for core analysis, including canonical pathways, functional analysis, and disease-related pathways.

Construction of the Protein–Protein Interaction (PPI) Network

To explore the potential interactions among the hub genes, the DE-IRGs were submitted to the Search Tool for the Retrieval of Interacting Genes (STRING) database. After removing the genes of degree < 2, the PPI network of DE-IRGs was visualized by Cytoscape software. Yellow nodes of this network represented hub genes.

Receiver Operating Characteristic (ROC) Curve and Logistic Regression Analysis

ROC curve analysis of the hub genes was performed to evaluate the sensitivity (true positive rate) and specificity (true negative rate) for distinguishing COPD patients. Thereafter, the area under the curve (AUC) of the ROC curve was calculated using the R package pROC. To further assess the efficacy of the hub genes in diagnosing COPD, we performed a logistic regression model in the GSE57148 dataset. The logistic regression model was constructed with the scikit-learn package of Python software.

Statistical Analysis

All statistical analyses were performed using the R version, and the data from different groups were compared by the Wilcoxon test. The correlations between immune cells and DE-IRGs were determined using Spearman correlation analysis. All P values of the present study were two-tailed, and P < 0.05 was considered statistically significant.

Results

The Profile of Immune Infiltration in COPD

We first revealed the profile of 22 infiltrating immune cell types in COPD, and subsequently, we analyzed the difference between COPD and normal samples using the CIBERSORT algorithm. Figure 1A shows the general distribution of the 22 immune cell types in each sample between the two groups. Compared with normal samples, lung tissue of COPD patients contained a greater number of resting NK cells, activated dendritic cells, and neutrophils. However, the fractions of follicular helper T cells and resting dendritic cells were relatively lower (Figure 1B, P < 0.05). Neutrophils were positively correlated with monocytes (Pearson correlation = 0.38, P < 0.01), while activated mast cells were negatively correlated with resting mast cells (Pearson correlation = −0.66, P < 0.01) (Figure 1C). Taken together, these results indicated that the immune response of COPD plays an important role in the progression of COPD.

Figure 1 The profile of immune infiltration in COPD. (A) Immune infiltrating cell ratio between normal samples and COPD patients. (B) Comparison of the difference between immune infiltrating cells in COPD and normal samples. (C) Correlation of each immune infiltrating cell. The value represents the correlation coefficient between immune cells (range −1 to 1) were shown in the upper right half. Immune cells with higher, lower, and same correlation levels were shown in red, purple, and white, respectively. Significant P-values for correlations between immune cells were shown in the lower left half, * for P < 0.05, ** for P < 0.01, ***P < 0.0001.

Identification of Differentially Expressed Immune-Related Genes (DE-IRGs)

To further identify the DE-IRGs between the normal and COPD samples, we comprehensively explored the gene expression microarray of the GSE57148 dataset. As shown in Figure 2A, there were 378 DEGs (|log2 FC| < 0.5 and FDR < 0.05) in lung tissues affected by COPD versus those in unaffected normal lung tissues in the GSE57148 dataset, of which 25 were upregulated and 353 were downregulated (Supplementary Table 1). To further identify the DEGs that were strongly associated with immunity, we intersected the 378 DEGs of the GSE57148 dataset with 1509 genes obtained from the ImmPort database. Thirty-eight DE-IRGs were obtained, including SOCS3, ISG15, TNFAIP3, GNLY, EIF2AK2, BST2, AHNAK, FGA, NFKBIZ, DMBT1, IL7R, PTGS2, CXCR1, CCL24, IL6, MIF, NFATC2, HRAS, DEFB1, SLIT2, CXCR2, AREG, INHBA, LIF, NAMPT, NPFF, RETN, S100A6, SCGB3A1, EGFR, EPOR, IGF2R, NR4A3, NRP2, TNFRSF10C, TNFRSF14, FCGR3B, and LAT (Figure 2B). In addition, the expression of 38 DE-IRGs is shown in a heatmap (Figure 2C).

Figure 2 Identification of differentially expressed immune-related genes. (A) Volcano map of differentially expressed genes. (B) Venn diagram of differentially expressed genes and immune-related genes. (C) Heat map of differentially expressed immune-related genes.

Functional Enrichment Analysis of DE-IRGs

To further demonstrate which biological functions and pathways these DE-IRGs are involved in, we performed GO and KEGG pathway analyses on DE-IRGs. As expected, the biological processes of DE-IRGs were primarily associated with B cell activation and the regulation of cytokine production involved in the immune response. For molecular functions and cellular components, these DE-IRGs were also involved in cytokine binding, cytokine receptor activity, cytokine receptor binding, cytokine activity, receptor ligand activity membrane microdomain, membrane raft, and secretory granule membrane (Figure 3A, Supplementary Table 2, P adjust < 0.05). In addition, the KEGG pathway analysis showed that the DE-IRGs were significantly enriched in natural killer cell-mediated cytotoxicity, Kaposi sarcoma-associated herpesvirus infection, C-type lectin receptor signaling pathway, PD-L1 expression and PD-1 checkpoint pathway in cancer, human cytomegalovirus infection, TNF signaling pathway, JAK-STAT signaling pathway, viral protein interaction with cytokine and cytokine receptor, and cytokine-cytokine receptor interaction (Figure 3B, Supplementary Table 3, P adjusted < 0.05). Surprisingly, DE-IRGs were associated with coronavirus disease COVID-19 in the development of COPD.

Figure 3 Functional enrichment analysis of differentially expressed immune-related genes. (A) Differential gene GO annotation visualization (Top 10) (bubble chart). (B) Differential gene KEGG annotation visualization (Top 10) (bubble chart). (C) Immune-related genes, differential gene enrichment, classical signaling pathways. (D) Immune-related differential gene enrichment of disease-related pathways.

The IPA results further revealed that these DE-IRGs activated several immune-related pathways, such as NF-κB signaling, the Th17 activation pathway, IL-8 signaling, and systemic lupus erythematosus in the B cell signaling pathway (Figure 3C, FDR < 0.01). Moreover, DE-IRGs were significantly enriched in cellular movement, immune cell trafficking, infectious diseases, connective tissue disorders, inflammatory disease, and skeletal and muscular disorders (Figure 3D, FDR < 0.01).

Identification of Hub Genes in COPD

A previous study revealed that T lymphocyte subset (Th1/Th2, Th17/Treg) imbalance is involved in the pathogenesis of COPD.¹⁷ To further identify the hub genes in COPD, we first analyzed the correlations between the 38 DE-IRGs and 5 differentially expressed immune cells (Figure 4A). Correlation analysis indicated that AHNAK and SLIT2 were negatively correlated with follicular helper T cells in the GSE57148 dataset of COPD (Figure 4B, Cor = −0.651, P = 4.86e-23; Cor = −0.667, P = 1.43e-24, respectively). However, TNFRRSF10C, CXCR1, CXCR2, and FCGR3B were positively correlated with neutrophils in COPD progression (Cor = 0.67, P = 8.21e-25; Cor = 0.853, P = 3.31e-52; Cor = 0.854, P = 2.05e-52; Cor = 0.833, P = 1.36e-47, respectively). The six DE-IRGs associated with immune cells were considered hub genes. To explore the interaction of these hub genes, a PPI network consisting of nodes and edges was constructed with the STRING database, which contained CXCR1, CXCR2, and FCGR3B (Figure 4C). These results indicated that there were complex interactions among the hub genes in the occurrence and development of COPD. In addition, the expression of six hub genes in COPD patients was greater than that in the normal group (Figure 4D, P < 0.05). Further, we verified the expression levels of the above six hub genes in the GSE76925 dataset (Supplementary Figure 1). Among them, the expression trends of AHNAK and FCGR3B were opposite to theirs in the GSE57148 dataset, which might be related to sample heterogeneity.

Figure 4 Identification of hub genes in COPD. (A) Correlation diagram of differentially expressed genes and differential immune cells. (B) Scatter plot of immunological correlation between differentially expressed genes and differential immune cells. (C) PPI network of hub genes. (D) Boxplots showed the expression of hub genes between normal samples and COPD patients.

The Diagnostic Value of Hub Genes for COPD Patients

To assess the diagnostic significance of these hub genes for COPD patients, we calculated the AUC value to evaluate the sensitivity and specificity of the hub genes (Figure 5A). The AUCs of the hub genes were 0.793, 0.688, 0.774, 0.661, 0.655, and 0.644, revealing their diagnostic value in COPD. To further assess the efficacy of the hub genes in diagnosing COPD, we performed a logistic regression model in the GSE57148 dataset (Figure 5B). The sensitivity and specificity of the ROC curve were 0.922 and 0.722, respectively. Importantly, the AUC of this model was 0.866, indicating that the prediction model based on the 6 hub genes has high accuracy and specificity for COPD diagnosis. Moreover, the diagnostic value of the six hub genes and the logistic regression model based on them in the GSE76925 dataset was displayed in Supplementary Figure 2.

Figure 5 The diagnostic value of hub genes for COPD patients. (A) ROC curves of hub genes. (B) ROC curve of regression model based on hub genes.

Construction of a Network with TF-Hub Gene Pairs

To further explore the correlations between TFs and hub genes, a regulatory network of TF-hub gene pairs was constructed by Cytoscape. As shown in Figure 6, a transcription regulatory network with 30 edges and 29 nodes was obtained for these hub genes. The results revealed that different hub genes were regulated by different TFs, and one TF could regulate multiple hub genes simultaneously.

Figure 6 The regulatory network of transcription factors and hub genes.

Discussion

Based on the CIBERSORT analysis, novel diagnostic biomarkers in COPD are essential for predicting the risk of developing the disease among people exposed to risk factors or for predicting the severity of the COPD course.

In this study, the CIBERSORT results revealed that the lung tissues of COPD patients contained a greater number of resting NK cells, activated dendritic cells, and neutrophils, and the fractions of follicular helper T cells and resting dendritic cells were relatively lower. Based on the theory of the immune mechanism of COPD in recent years, a variety of new immunomodulatory drugs can promote damaged airway repair or reconstruction of immunity, thereby reducing the frequency and severity of COPD attacks and improving health status and exercise tolerance. Therefore, this study found differential immune cells in COPD samples, which can provide new ideas for the treatment of COPD.

Six DE-IRGs related to immune cells were identified as diagnostic genes for COPD, including AHNAK, CXCR1, SLIT2, CXCR2, TNFRSF10C and FCGR3B. AHNAK is an exceptionally giant protein (700 kDa) that was initially identified in human neuroblastomas and skin epithelial cells¹⁸ and regulates cytoskeletal structure formation, muscle regeneration, calcium homeostasis, and signaling.¹⁹ AHNAK functions as a tumor suppressor in murine lungs by suppressing alveolar epithelial cell proliferation and modulating the lung microenvironment.^20,21 In recent years, more attention has been given to the role of the AHNAK family in tumor progression.²² Upregulation of AHNAK was significantly associated with poor prognosis of laryngeal carcinoma, mesothelioma, and pancreatic ductal carcinoma.^23,24 AHNAK is a tumor suppressor in breast cancer due to its ability to activate the TGFb signaling pathway.²¹ Using the powerful genome-wide linkage scan in a Dutch genetic isolate, Ivana et al confirmed the implication of the 15q25 region in COPD and identified regions at chromosomes 5 and 11. Within the region on chromosome 11, four deleterious rare variants in AHNAK, PLCB3, SLC22A11 and MTL5 were identified that were shared among most of the affected family members.²⁵

SLIT2 is a secreted glycoprotein of the SLIT family that is involved in the epithelial-mesenchymal transition process that permits cancer cells to acquire migratory, invasive, and stem-like properties.^26–28 A previous study reported that SLIT2 inhibited CXCL12/CXCR4-induced chemotaxis, invasion, adhesion of breast cancer cells and the secretion of MMP-2 and MMP-9.²⁹ SLIT2 was downregulated in cigarette smoke-exposed cells and lung adenocarcinoma, and increased expression of SLIT2 was associated with a better outcome in lung adenocarcinoma patients.³⁰ We also found that AHNAK and SLIT2 were negatively correlated with follicular helper T cells in the GSE57148 dataset of COPD.

In our study, we found that TNFRRSF10C, CXCR1, CXCR2, and FCGR3B were positively correlated with neutrophils in COPD progression. CXCR1 and CXCR2 can promote neutrophil migration, and their chemokine receptors are G protein-coupled receptors that bind to human IL-8 family chemokines (CXCL1, CXCL2, CXCL3, CXCL5, CXCL6, CXCL7 and CXCL8).³¹ CXCR1 and CXCR2 were observed on NK cells and lung structural cells, such as fibroblasts and epithelial and endothelial cells. An increasing number of studies have found that CXCR1 and CXCR2 play a role in the pathogenesis of certain pulmonary diseases, such as COPD, asthma, and lung fibrosis.^32,33 Immunostaining for CXCR1 was observed in inflammatory cells, endothelial cells and, to a lesser extent, in fibroblasts within the bronchial submucosa.³⁴ CXCR1 activates neutrophils and mediates neutrophil functions, including chemotaxis, intracellular calcium changes, phospholipase D activation, and respiratory burst, through the binding of specific ligands such as IL-8.^35–37 CXCR1 expression in circulating neutrophil cells was significantly higher in COPD patients than in healthy subjects because the overexpression of CD-11b and CXCR1 in circulating neutrophils may be associated with the development of airflow limitation in COPD patients.³⁸ The researchers found that the amount of smoking correlated only with CXCR1 expression in subjects who had a history of smoking (patients with COPD and smokers). Long-term exposure to smoking may play some role in the regulation of CXCR1.³⁸ CXCR2, a G-protein-coupled receptor, binds human CXC chemokine ligands. CXCR2 is upregulated in exacerbations of COPD, where its expression colocalizes with the accumulation of airway mucosal neutrophils.³⁹ Therefore, a study showed that the CXCR2 antagonist was effective in COPD.⁴⁰ Blood neutrophil counts were associated with exacerbation frequency and mortality in COPD.⁴¹ IL-8 exerts its potential effects on neutrophils by binding with high affinity to the chemokine receptors CXCR1 and CXCR2 on its cell surface.⁴² Keir et al demonstrated that there was a subset of individuals with COPD in whom neutrophil activation was independent of CXCR2.⁴³ Another study found that specific CXC chemokines resulted in enhanced migration of PBMCs/monocytes from patients with COPD, which was mediated via CXCR1 and CXCR2.⁴⁴ CXCR2 was significantly increased in patients with severe/very severe COPD compared with patients with mild/moderate COPD and control healthy smokers.³⁴

TNFRSF10C, a member of the TNF family, is a decoy receptor for TNF-related apoptosis-inducing ligand (TRAIL).⁴⁵ Hypermethylation of TNFRSF10C was associated with non-small cell lung cancer (NSCLC).⁴⁶ TNFRSF10C hypermethylation was significantly associated with the risk of colorectal cancer.⁴⁷ We know that is associated with some immune system diseases. For example, FCGR3B is associated with rheumatoid arthritis, and this gene may play a role in the physiopathology of this disease.⁴⁸ A low copy number of FCGR3B is a risk factor for lupus nephritis in a Chinese population.⁴⁹ However, TNFRSF10C and FCGR3B have not been reported in COPD. This is the first time that we have reported their diagnostic value in COPD.

Functional enrichment analysis of 38 DE-IRGs identified enrichment of immune-related pathways. Enriched KEGG pathways included PD-L1 expression and the PD-1 checkpoint pathway in cancer and other disease pathways, which were all related to COPD. In our future studies, we will perform research on the immunotherapy of COPD according to these immune-related pathways and the six identified diagnostic genes. These DE-IRGs were also associated with COVID-19 in the development of COPD, which was consistent with a previous study. COPD patients have an increased risk of severe pneumonia and poor outcomes when they develop COVID-19.⁵⁰ This provides a new idea for the study of COPD and COVID-19.

Conclusion

AHNAK, SLIT2, TNFRRSF10C, CXCR1, CXCR2, and FCGR3B may have a strong influence on COPD. These hub genes may serve as control targets for COPD and further study of these genes may lead to new therapies.