Immune Characteristics Analysis and Transcriptional Regulation Prediction Based on Gene Signatures of Chronic Obstructive Pulmonary Disease

Highlights

1. Key module-related DEGs in COPD were found to involve in interleukin-related pathways.

2. An optimized gene set, including seven-gene signatures, was identified using LASSO analysis.

3. The seven-gene signature was significantly correlated with 18 types of immune cells in COPD.

4. The upstream miRNAs and TFs of gene signatures were predicted.

Introduction

Chronic obstructive pulmonary disease (COPD) is a complex heterogeneous disease characterized by persistent respiratory symptoms and airflow limitation because of respiratory and/or alveolar abnormalities.^1,2 A forced expiratory volume in 1s (FEV1)/forced vital capacity (FVC) less than 0.70 after using bronchodilators is a necessary condition for the diagnosis of COPD.³ It is estimated to cause more than three million deaths worldwide annually, affecting about one in 10 adults aged over 45.^4,5 Although smoking is the leading cause of COPD, the unsatisfactory prevalence of COPD in non-smokers suggests that chronic inflammation in patients may be caused by different cellular and pathophysiological changes in different genetic backgrounds.⁶ Some studies believe that genetic risk factors contribute to the development of airflow restriction, and the huge variation in lung inflammation between individuals is also genetically determined.⁷ Therefore, effective biomarkers are needed to predict risk and guide personalized treatment of patients with COPD.

Biomarkers are of great significance in the diagnosis, prognosis, and treatment of diseases. In recent years, many clinically relevant COPD biomarkers, including plasma fibrinogen, CRP, IL-6, IL-8, total bilirubin, SAA, SP-D, CCSP-16, MMP-8, and MMP-9, have been developed using the candidate gene approach.^8–12 However, the influence of a single biomarker is often limited in predicting the disease. It has been reported that multiple biomarkers are more powerful than a single biomarker, and their combination with clinical variables can better predict disease progression.¹³ Based on shared data from high-throughput sequencing and microarray technology, it is possible to observe changes in the expression of numerous genes in different samples at one time.¹⁴ Moreover, bioinformatics methods that integrate multiple datasets can help identify reliable and repeatable genetic signatures associated with disease development.¹⁵ Therefore, data from four microarrays were used to analyze differentially expressed genes (DEGs) in this study. Then, a weighted gene co-expression network analysis (WGCNA) and the least absolute shrinkage and selection operator (LASSO) regression analysis were performed to obtain the optimized gene set from the trait-related module, which was also shown to have a great predictive ability among COPD patients.

Additionally, inflammation in the lungs of patients with COPD is associated with abnormal immune responses, and both the innate and adaptive immune systems are involved, linked by the activation of dendritic cells.¹⁶ Furthermore, an increase in the number of macrophages, neutrophils, and T and B lymphocytes has also been observed in the COPD process, and they are recruited in the immune cycle. They can secrete proinflammatory factors, such as cytokines, chemokines, growth factors, and lipid mediators.^17,18 Considering the key roles of these immune cells, we attempted to identify the changes in the immune microenvironment defined by the gene signatures selected above and explain the molecular regulation mechanisms of upstream transcription factors (TFs) and miRNAs. A flowchart of this study is presented in Figure S1.

Methods

Data Acquisition and Process

All data were downloaded from the Gene Expression Omnibus database (GEO, http://www.ncbi.nlm.nih.gov/geo/). Using COPD, lung, and smokers as keywords, the relevant data of humans in the last five years were searched, and a total of four microarrays were obtained, including GSE37768, GSE106986, GSE103174, and GSE76925. Among them, GSE37768 included 18 lung tissue samples from COPD patients and nine healthy smokers and non-smokers, which were detected by GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array; GSE106986 included 14 lung tissue samples from COPD patients and five healthy non-smoker samples, which were detected by GPL13497 Agilent-026652 Whole Human Genome Microarray 4×44K v2; GSE103174 included 37 lung tissue samples from COPD patients and 10 healthy non-smoker samples which were detected by GPL13667 [HG-U219] Affymetrix Human Genome U219 Array; GSE76925 included 111 COPD samples and 40 smoker control samples which were detected by GPL10558 Illumina HumanHT-12 V4.0 expression beadchip. Among them, GSE76925 was used for the external validation. The detailed clinical information of the samples from these four datasets is summarized in Table S1.

Data Pre-Processing

The standardized probe expression value matrices of the above datasets were downloaded from GEO, and annotation files were downloaded from the corresponding platforms. Probes that did not match the gene symbol were removed, and when different probes were mapped to the same gene, the mean value of the probes was taken. Based on GSE37768, GSE106986, and GSE103174, the ComBat function in R sva (version 3.34.0, http://www.bioconductor.org/packages/release/bioc/html/sva.html)^19,20 was used to remove the batch effect, and a combined dataset was obtained for further analysis.

DEGs Analysis

To analyze the expression differences between COPD samples and controls in the combined dataset, the Limma package (version 3.10.3, http://www.bioconductor.org/packages/release/bioc/html/limma.html)²¹ was used. The p-value obtained after the t-test was corrected by multiple tests using the Benjamini & Hochberg (BH) method, and the results were then evaluated in terms of difference multiples and significance. Adjusted p-value < 0.05, and |log fold change (FC)| > 0.585 were set as standards to select DEGs.

WGCNA

Genes with an absolute deviation of median expression in the top 2000 were selected for the WGCNA algorithm (version 1.61, https://cran.r-project.org/web/packages/WGCNA/).²² WGCNA defined the correlation matrix of gene co-expression and the power exponential linkage function of the gene network, and the power value was captured when the square value of the correlation coefficient reached 0.9 for the first time. Then, a hierarchical clustering tree was constructed by calculating the dissimilitude degree of different nodes. The module was screened with a minimum number of genes in each gene network set to 30 and a pruning cut height of 0.3. Finally, the trait-related module was identified by analyzing the correlation between the module and disease phenotype with standards of p < 0.05 and |r| > 0.3. By considering the intersection of DEGs and genes in trait-related modules, key module-related DEGs were screened out using the R package VennDiagram (version 1.6.20, https://CRAN.R-project.org/package=VennDiagram).²³

Enrichment Analysis

Based on the key module-related DEGs obtained above, over-representation analysis in the functional profiling module of gprofiler online tools (https://biit.cs.ut.ee/gprofiler/convert)²⁴ was performed for gene ontology (GO) functions²⁵ and pathways²⁶ enrichment analyses. GO functions included biological process (BP), cellular component (CC), and molecular function (MF), whereas pathways were enriched in the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome databases. BH adjusted p < 0.05 was considered statistically significant.

Screening and Validation of Gene Signatures

The key module-related DEGs were then incorporated into the LASSO regression model using the glmnet package (version 4.0–2, https://cran.r-project.org/web/packages/glmnet/index.html),²⁷ and the gene signatures were obtained with 20-fold cross-validation. Gene signatures were further verified on the combined datasets and GSE76925 by constructing a support vector machine (SVM) model using the e1071 SVM package (version 1.7–4, https://CRAN.R-project.org/package=e1071).²⁸ Receiver operator characteristic (ROC) curves were also created to verify the accuracy of the model.

Immune Cell Infiltration Analysis

The CIBERSORT deconvolution algorithm (https://cibersort.stanford.edu/index.php)²⁹ was used to estimate the infiltration abundance of 22 types of immune cells in each sample of the combined datasets, with the parameters set as perm = 100 and QN = F. Then, the Wilcoxon test was used to analyze the differences in the abundance of immune cell infiltration between the COPD and control groups, and a violin plot was created using the vioplot in R package (version 0.3.5, https://github.com/TomKellyGenetics/vioplot).³⁰ Finally, the Pearson correlation coefficient of gene signatures and immune cell infiltration abundance was calculated, and relation pairs (p < 0.05) were selected.

Predicting Upstream TFs and miRNA of Gene Signatures

The relationships between miRNA and gene signatures were predicted using miRWalk 3.0 (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk3.html).³¹ The relationship pairs appearing in both miRDB (http://www.mirdb.org/miRDB/download.html)³² and miRTarBase (http://mirtarbase.cuhk.edu.cn/php/index.php)³³ databases with scores over 0.9 were selected to construct a miRNA regulatory network. The TH-mRNA relations were predicted in the TRRUST database (https://www.grnpedia.org/trrust/),³⁴ and the TF-mRNA and miRNA-mRNA relation pairs were integrated into a TF-mRNA-miRNA network.

Results

Screening of DEGs Between COPD and Healthy Controls

To remove the batch effect between GSE37768, GSE106986, and GSE103174, a principal component analysis was performed, and a combined dataset was obtained with thresholds of |logFC| > 0.585 and adjusted p < 0.05. A total of 127 DEGs between COPD and control samples in the combined dataset were screened out, as shown in Figure 1A. Detailed information on these DEGs is presented in Table S2. Among them, 104 DEGs were up-regulated, and 23 DEGs were down-regulated (Figure 1B).

Figure 1 Screening of DEGs between COPD and healthy controls. (A) The heatmap shows the DEGs between COPD and control samples from the combined dataset of GSE37768, GSE106986, and GSE103174, with thresholds of |logFC| > 0.585 and adjusted p < 0.05. (B) Volcano plot showing up-regulated (red) and down-regulated (blue) DEGs. The top five up-regulated and down-regulated DEGs ranked by |logFC| are labeled in red and blue, respectively. Abbreviations: COPD, chronic obstructive pulmonary disease; FC, fold change.

WGCNA to Identify Trait-Related Modules and Key Genes

The WGCNA algorithm was then used to select trait-related modules and key genes. To satisfy the premise of a scale-free network distribution, the value of the power (weight parameter) of the adjacency matrix was explored. The selection range of network construction parameters was first set, then the scale-free distribution topology matrix was calculated, and the graph was created as shown in Figure 2A. When the square value of the correlation coefficient reaches 0.9 for the first time, the value of power is 6. A systematic cluster tree was obtained by calculating the coefficient of dissimilarity between genes. Then, the dynamic mixing cutting method was used to divide the modules, and finally, eight disease-related modules were obtained and labeled with different colors (Figure 2B). The genes involved in each module are listed in Table S3. Thereafter, the correlation between modules and disease traits was evaluated by calculating the correlation coefficient between the gene significance of each module and the diseased trait. The results (Figure 2C) suggested that module turquoise and module yellow had the highest correlation with disease traits, and these two modules were also identified as trait-related modules. By considering the intersection of genes in trait-related modules and DEGs obtained above, a total of 83 common genes were identified as key module-related DEGs (Figure 2D, Table S4). Among them, 16 DEGs in module yellow were down-regulated, and 67 DEGs in module turquoise were up-regulated.

Figure 2 WGCNA algorithm was conducted to identify trait-related modules and genes. (A) Selection of power, weight parameter of the adjacency matrix. The x-axis represents the soft threshold (power), whereas the y-axis represents the square of the correlation coefficient between log2k and log2p(k) in the network. The red line represents the standard line, where the square value of the correlation coefficient reaches 0.9. (B) Module partition tree. The upper diagram represents the systematic cluster tree, and the different colors in the lower dynamic tree cut represent different gene modules. (C) Correlation between each module and disease traits. Statistical significance was set at p < 0.05. (D) The Venn plot shows the key module-related DEGs by taking the intersection of genes in trait modules and DEGs in COPD. Abbreviation: DEG, differential expressed gene.

Functions and Pathways Enrichment Analyses

Based on the key module-related DEGs, enrichment analyses were performed to investigate their functions and pathways they may be involved in. As a result, 626 GO-BP, 13 GO-MF, 22 KEGG pathways, and 12 Reactome pathways were enriched. The top five results of each term ranking by p-value are shown in Figure 3A and B. It can be observed that these key module-related DEGs were mainly enriched in GO-BP of the response of cytokine and cellular response to cytokine stimulus; GO-MF of cytokine and chemokine activities; KEGG pathways of IL-17 and TNF signaling pathways; and Reactome pathways of IL-4, IL-13, and IL-10 signaling pathways. The genes involved in these key functions and pathways are listed in Table S5.

Figure 3 Functions and pathways enrichment analyses. (A) Top5 GO functions enriched in BP and MF. (B) Top5 KEGG pathways and Reactome pathways were enriched. The results of each term were sorted using p-values. Abbreviations: GO, gene ontology; MF, molecular function; CC, cellular component; BP, biological process; KEGG, Kyoto Encyclopedia of Genes and Genomes; REAC, Reactome; IL, interleukin; TNF, tumor necrosis factor.

Identification of Gene Signatures

LASSO regression analysis was further performed to select the optimized gene set, and the parameters of the LASSO model are shown in Figure 4A. As a result, MTHFD2, KANK3, GFPT2, PHLDA1, HS3ST2, FGG, and RPS4Y1 were included in the optimized gene set when the lambda value was close to lambda.1se, and were thus identified as gene signatures. The expression of these gene signatures was then validated in the combined dataset (Figure 4B), and the results suggested that they were all significantly differentially expressed between COPD and healthy samples. ROC curves were also created to evaluate the efficacy of the seven-gene-signature in predicting COPD, as shown in Figure 4C and D. The results suggested that the area under the curve (AUC) created for the combined dataset was 0.795, whereas the AUCs of GSE37768, GSE106986, GSE103174, and GSE76925 were all over 0.700, thereby indicating a good predictive ability of gene signatures in COPD.

Figure 4 Identification of gene signatures using LASSO regression analysis. (A) The optimized gene set with seven DEGs was selected as the gene signature when the lambda value was close to lambda.1se. (B) Validation of gene signatures in the combined dataset. ****p < 0.0001. (C) ROC curve created for the combined dataset. (D) The ROC curves of GSE37768, GSE106986, GSE103174, and GSE76925. Abbreviations: ROC, receiver operator characteristic; SVM, support vector machine; AUC, area under the curve.

Correlation Analysis of Gene Signatures and Immune Cell Infiltration Abundance

As described in the Methods section, the infiltration abundances of 22 types of immune cells were calculated for each sample. The differences in infiltration abundance between the COPD and control samples were analyzed, as shown in Figure 5A. The results suggested that fractions of plasma cells, CD8 T cells, resting NK cells, activated dendritic cells, and activated mast cells were significantly different between the two groups (p < 0.05). Thereafter, the correlation coefficients between gene signatures and immune cells were calculated, and relationships with p < 0.05 were considered statistically significant. As shown in Figure 5B and Table 1, FGG, GFPT2, HS3ST2, KANK3, MTHFD2, PHLDA1, and RPS4Y1 were significantly correlated with six, eight, eight, six, eleven, ten, and four types of immune cells. In these 53 gene-cell relationships, plasma cells were significantly related to MTHFD2, GFPT2, PHLDA1, HS3ST2, and RPS4Y1; CD8 T cells were significantly related to six gene signatures except PHLDA1; resting NK cells were significantly related to six gene signatures except RPS4Y1; activated dendritic cells were significantly related to FGG and PHLDA1; and activated mast cells were significantly related to five gene signatures except KANK3 and RPS4Y1.

Table 1 Correlations Between Gene Signatures and Immune Cells

Figure 5 Correlation analysis of gene signatures and immune cells. (A) The violin plot shows differences in immune cell fractions between COPD and control samples. Yellow and blue bars represent the cases and controls, respectively. *0.01<p < 0.05; **0.001<p < 0.01; ***p < 0.001. (B) Correlation analysis of seven-gene signatures and 21 types of immune cells. Statistical significance was set at p < 0.05.

Construction of a miRNA-mRNA-TF Network

In the miRWalk database, a total of nine upstream miRNAs were predicted to target with three gene signatures, including GFPT2, MTHFD2, and PHLDA1, and nine miRNA-mRNA relation pairs were obtained. Furthermore, two TFs were predicted to target PHLDA1 and FGG, and three mRNA-TF relation pairs were obtained. By integrating miRNA-mRNA and mRNA-TF relations, a miRNA-mRNA-TF network comprising 16 nodes and 12 regulatory axes was established, as shown in Figure 6. Notably, all mRNAs in this network were up-regulated DEGs in the turquoise module. Among them, PHLDA1 was identified as the hub gene that interacts with the miRNA (miR4316) and the TFs (FLI1 and EWSR1).

Figure 6 Construction of a miRNA-mRNA-TF network by integrating miRNA-mRNA and mRNA-TF relations. The turquoise, red, and purple dots represent mRNA, miRNA, and TFs, respectively.

Discussion

In this study, we first screened out 127 DEGs between COPD and control samples from the combined dataset of GSE37768, GSE106986, and GSE103174. By considering the intersection of 127 DEGs and genes in trait-related modules, 83 key module-related DEGs were identified, which were mainly enriched in interleukin-related pathways. Then, an optimized gene set, including MTHFD2, KANK3, GFPT2, PHLDA1, HS3ST2, FGG, and RPS4Y1, was obtained using the LASSO regression analysis, and it showed great potential in predicting COPD risks in the combined dataset with an AUC of 0.795.

The effect of the seven-gene-signature on COPD has been limited. Among them, MTHFD2 is associated with the expression of cell cycle-related genes in non-small cell lung cancer (NSCLC) and has clinical significance in the diagnosis, pathological stage, and prognosis of the disease.^35,36 During tumor development, MTHFD2 interacts with a series of nucleoproteins involved in RNA metabolism and translation.³⁷ In predicting upstream miRNAs in this study, we found that MTHFD2 may have potential regulatory axes with miR-9-5p, miR-6880-5p, miR-1910-5p, and miR-3165, of which the expression relationship between MTHFD2 and miR-9 has been confirmed in breast cancer.³⁸ Therefore, we hypothesized that the miRNA-MTHFD2 regulatory axis might affect cell carcinogenesis in COPD patients by regulating the expression of cell cycle-related genes. PHLDA1 plays an important role in cell proliferation and apoptosis and is involved in lung contusion and subsequent inflammatory responses.³⁹ In this study, it was found that PHLDA1 was significantly highly expressed in patients with COPD, and the expression of PHLDA1 was also shown to be significantly correlated with the prognosis of NSCLC,⁴⁰ indicating a significant contribution of PHLDA1 in the process of NSCLC from COPD. Although we have proposed that candidate genes may be involved in the progression of COPD to NSCLC, whether these genes affect the severity of COPD remains unclear. This study also predicted the potential transcriptional regulation of PHLDA1 with FLI1 and EWSR1, whereas Boro et al supported our prediction and stated that the expression of PHLDA1 might be inhibited by EWS/FLI1.⁴¹

By comparing the differences in infiltration abundance of immune cells between COPD and controls, we found that the fractions of plasma cells, CD8 T cells, resting NK cells, activated dendritic cells, and activated mast cells were significantly different between two groups. Studies have confirmed the interactions between mast cells, airway structure cells, and other inflammatory cells in COPD development.⁴² In addition, perivascular mast cell density is positively correlated with airway angiogenesis in COPD, which contributes to airway remodeling.⁴³ Furthermore, mast cells were found to accumulate abnormally in COPD and may increase the expression of IL-17.⁴⁴ Interestingly, the key module-related genes obtained in this study were also enriched in IL-17 related signaling pathways. The infiltration of dendritic cells was significantly increased in patients with COPD in the current study. Van et al agreed with this finding and stated that in patients with COPD, especially in severe disease, dendritic cells were activated and increased in number.⁴⁵ Van et al also proposed that the number of Langerhans-type dendritic cells increases with COPD severity.⁴⁶ It has been reported that immature dendritic cells in the respiratory tract promote Th2 differentiation.⁴⁷ Moreover, the epithelial cells of asthmatic patients can attract Th2 cells into the airway by releasing chemokines through dendritic cells.⁴⁸ By studying the relationship between gene signatures and dendritic cells, we found that PHLDA1 and FGG were positively correlated with dendritic cells, and PHLDA1 and FGG were highly expressed in COPD patients. These results suggest that PHLDA1 and FGG are up-regulated in COPD mediated by TFs and may be involved in the immune regulation mechanism associated with dendritic cells, thereby influencing the progression of COPD.

It is undeniable that the integration of mRNA and non-coding RNA research techniques provides a more comprehensive biological perspective for exploring the pathogenesis of diseases. However, the relationship between gene signatures and miRNAs was only investigated at a predictive level in this study. The lack of exploration of the expression correlation and regulation mechanism between them is one of the limitations of this study. Additionally, the COPD severity stratification is unavailable in this study because of the lack of clinical features. Therefore, in the follow-up study, we will collect clinical solid tumor samples and further include information on disease severity to explore the role of genetic markers in COPD progression. Further functional experiments will be conducted to confirm the regulatory relationship between miRNA-mRNA and TF-mRNA regulatory axes.

Conclusion

To conclude, this study proposed seven-gene signatures of COPD, including MTHFD2, KANK3, GFPT2, PHLDA1, HS3ST2, FGG, and RPS4Y1, based on difference analysis and trait module analysis. These genes can be used as potential diagnostic markers to predict the onset risk of COPD. The seven-gene signature was enriched in interleukin-related signaling pathways and was also associated with immune cells. Additionally, the miRNAs and TFs of gene signatures were analyzed to help understand upstream regulatory mechanisms. This study provides valuable biomarkers for patients with COPD to customize personalized treatments.