Identification of Dysregulated Mechanisms and Candidate Gene Markers in Chronic Obstructive Pulmonary Disease

Introduction

According to estimates from the Global Burden of Disease (GBD) study, 523 million around the world had cardiovascular disease in 2019.¹ Chronic obstructive pulmonary disease (COPD) is a major cause of global morbidity and mortality related to cardiovascular disease. COPD is caused by persistent, often progressive airflow limitation in the lungs, and it includes the conditions of chronic bronchitis and emphysema.² More than 3 million people die from COPD every year worldwide,³ and COPD is expected to become the third leading cause of death globally by 2030.⁴ Despite the worldwide prevalence of this disease, it remains largely underdiagnosed and undertreated.

Many factors have been associated with COPD, including systemic and local inflammation, air pollution and smoking.⁵ However, the exact mechanism of COPD remains unclear. Clinically, this disorder is characterized by cough and chronic dyspnea resulting from airway stenosis. Accumulating evidence indicates that COPD is a systemic disease, and that its pathological manifestations are not limited to lung inflammation and airway remodeling⁶ In fact, COPD has a profound effect on cardiac function and gas exchange, and COPD patients are at 2–5 times greater risk of cardiovascular disease than the general population.⁷ Cardiovascular disease in any stage of COPD greatly increases the risk of death and hospitalization.⁸ Approximately 22–40% of COPD patients experience at least one moderate or severe exacerbation each year, and 9–16% experience more than one exacerbation per year.⁹

Dyspnea is the main symptom of COPD and the most frequent reason why COPD patients seek medical attention.¹⁰ At the onset of COPD, timely intervention to relieve symptoms and exacerbations can prevent an acute decline in lung function and progression of the condition to severe emphysema.¹¹ Current treatments of COPD involve long-term inhalation therapy with bronchodilators, corticosteroids, or a combination of these agents.¹² Combination therapy may be more effective than monotherapies for relieving COPD exacerbations, but the evidence is inconclusive.¹³ A challenge to effective treatment is that the nine approved drug classes for COPD maintenance therapy treat only the symptoms, rather than the underlying inflammation or progression.¹⁴

Elucidating the pathogenic pathways in COPD may help identify new treatments and strategies to prevent exacerbations. Toward that end, we applied bioinformatics to public databases to explore molecular mechanisms and feature genes potentially associated with COPD. We explored the potential functions of the feature genes in the disease, with a focus on immune cell infiltration as well as gene regulation by microRNAs (miRNAs) and DNA methylation. Understanding these regulatory mechanisms and identifying potential marker genes may lead to the development of new therapeutic strategies.

Materials and Methods

Data Collection

All data in this study were obtained from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). The GSE47460 dataset¹⁵ included mRNA expression profiles from whole-lung homogenates of 144 COPD patients and 91 controls, obtained using the GPL14550 platform. The GSE151052 dataset¹⁶ included mRNA expression profiles from lung tissues of 77 COPD patients and 40 controls, obtained using the GPL17556 platform. The GSE38974 dataset¹⁷ included mRNA expression profiles from lung tissues of 23 COPD patients and 9 controls, obtained using the GPL4133 platform; as well as miRNA expression profiles from 19 COPD patients and 8 controls, obtained using the GPL7723 platform. The GSE76925 dataset¹⁸ included mRNA expression profiles from lung tissues of 111 COPD patients and 40 controls, obtained using the GPL10558 platform. The GSE63704 dataset¹⁹ included methylation profiles from lung tissues of 86 COPD patients and 26 controls, obtained using the GPL13534 platform.

Analysis of Expression and Methylation Differences Between COPD Patients and Controls

The limma package in R²⁰ was used to search for differences in mRNA levels between COPD patients and controls in the GSE47460 and GSE151052 datasets, as well as differences in miRNA levels between COPD patients and controls in the GSE38974 dataset. The cutoff value for differential gene expression was P < 0.05. DEmRs and DEmiRs were defined as differentially expressed mRNAs or miRNAs, respectively. DEmRs that were differentially expressed (up- or downregulated) in the same direction in both the GSE47460 and GSE151052 datasets were defined as common genes.

The cAMP package in R was used to identify differences in gene methylation between COPD patients and controls in the GSE63704 dataset. Differentially methylated probes (DMPs) were defined as methylation sites that passed the significance threshold of P < 0.05. Only DMPs with ∆Beta values that varied between patients and controls in the opposite direction as common genes were retained in the analysis.

Analysis of Functional Enrichment

Common genes were analyzed for enrichment in Gene Ontology (GO) biological processes and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways using the clusterProfiler package in R.²¹ Enrichment was considered significant at P < 0.05. Activated or inhibited enrichment was assessed using gene set variation analysis (GSVA) in R.²² Gene set enrichment analysis (GSEA) was also carried out using the fgsea package in R.

Construction of a Gene Coexpression Network

Common genes were analyzed using multiscale embedded gene coexpression network analysis (MEGENA) in R.²³ This analysis generates embedded, multi-scale networks to uncover biologically meaningful genes, which are then assigned to modules of coexpressed genes.

Regression Based on the Least Absolute Shrinkage and Selection Operator (LASSO)

Genes in modules 2, 3, and 6 were assessed for their ability to differentiate COPD patients and controls based on the area under the receiver operating characteristic curve (AUC), as calculated using the pROC package in R.²⁴ DEmRs that were present in both GSE47460 and GSE151052 datasets and that gave AUC > 0.7 were used to build a binomial LASSO regression model in the glmnet package in R.²⁵ LASSO regression was performed to shrink the regression coefficients towards zero as λ increased. We then optimized λ in order to identify feature genes.

Immune Cell Infiltration

Levels of infiltration by different types of immune cells in lung tissues of COPD were evaluated using single‐sample GSEA as coded in the GSVA package in R. The analysis was performed on data from datasets GSE151052, GSE38974, GSE47460, and GSE76925. Differences in infiltration levels between COPD patients and controls were identified using the limma package. We also evaluated potential correlations between feature gene expression and levels of immune cell infiltration using Pearson correlation analysis. Results associated with P < 0.05 were considered significant. Proportions of immune cells in lung tissues of COPD patients were analyzed using CIBERSORT (https://cibersort.stanford.edu/), after excluding immune cell types associated with “0”.

Target Prediction

Downstream genes possibly regulated by DEmiRs were predicted using Targetscan (http://www.targetscan.org/vert_72/).

Statistical Analysis

All statistical analyses were performed using R (version 3.9.1) and related packages within the R environment. Results associated with P < 0.05 were considered statistically significant.

Results

DEmRs in COPD

The design of this study is shown in Figure 1. We identified 6544 DEmRs between COPD patients and controls in the GSE47460 dataset (Figure 2A, Table S1) and 10,693 DEmRs in the GSE151052 dataset (Figure 2B, Table S2), based on a definition of differential expression as fold change > 0. A total of 487 DEmRs were upregulated in both datasets (Figure 2C), while 863 were downregulated in both (Figure 2D). These DEmRs were defined as common genes, which may be linked to COPD disease.

Figure 1 Flowchart of this study. Abbreviations: AUC, area under the receiver operating characteristic curve; COPD, chronic obstructive pulmonary disease; GSEA, gene set enrichment analysis; LASSO, least absolute shrinkage and selection operator.

Figure 2 Differentially expressed mRNAs between COPD patients and controls. (A) Differentially expressed mRNAs between COPD patients and controls in the GSE47460 dataset. Red dots are upregulated genes; green dots, downregulated genes. (B) Differentially expressed mRNAs between COPD patients and controls in the GSE151052 dataset. Venn diagram of (C) upregulated or (D) downregulated genes common to both datasets.

Biological Functions of Common Genes

To predict the biological functions in which common genes may be involved, we performed enrichment analysis. Enrichment analysis implicated the common genes in the GO biological processes of T cell activation, regulation of leukocyte activation, and lymphocyte differentiation (Figure 3A). The common genes were also implicated in the KEGG signaling pathways of leukocyte transendothelial migration, primary immunodeficiency, and asthma (Figure 3B). GSEA of common genes suggested that COPD was associated with activation of the intestinal immune network involved in IgA production, with activation of asthma pathways, and with inhibition of axon guidance and tight junctions (Figure 3C). These results suggest that common genes are involved in biological processes related to immune and inflammatory responses.

Figure 3 Potential functions of common genes in COPD. (A) Gene Ontology biological processes enriched in common genes. (B) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched in common genes. (C) The fgsea results of activated or inhibited KEGG pathways in COPD. NES, normalized enrichment score; pval, P value; padj, adjusted P value.

Construction of a Gene Coexpression Network

To further identify aberrantly expressed genes with important roles in COPD, we performed MEGENA on common genes (Figure 4A), which defined 186 modules containing a total of 1345 genes. The three largest modules were C1_ 6, with 270 genes; C1_3, 269 genes; and C1_2, 259 genes (Figure 4B). Genes in these three modules that differentiated COPD patients from controls with an AUC > 0.7 in both the GSE47460 and GSE151052 datasets (Figure 4C) were defined as candidate genes. We obtained 22 candidate genes (Figure 4D), which may be useful in diagnosing COPD.

Figure 4 MEGENA to identify COPD candidate genes based on coexpression of common genes. (A) Global MEGENA network of common genes. Nodes represent different modules. The larger the node, the greater the number of genes in the module. (B) Child modules with the largest number of genes in the MEGENA network. Different colors represent different child modules, and triangles represent hub genes of modules. (C) Areas under the receiver operating characteristic curve (AUC) of genes in modules C1_6, C1_3, and C1_2 in the GSE47460 and GSE151052 datasets. (D) Heatmap of candidate genes, where red indicates upregulated and green indicates downregulated.

Identification of Feature Genes in COPD

We used LASSO regression to select which of the 22 candidate genes were most likely to be important in COPD, based on the GSE47460 dataset. Based on an optimized λ of 0.02161155, we obtained 11 feature genes with non-zero coefficients (Figure 5A and B): AQP1, CLEC4A, GPC4, GPR39, PDCL3, PGM1, PRRG4, ROR1, RS1, SEC14L4, and USP13. After randomly dividing up the COPD patient samples in the GSE47460 dataset into a training set (75%) and validation set (25%), we evaluated the diagnostic performance of the feature genes against the training set based on AUC (Figure 5C). The AUC for all feature genes together was 0.885 against the training set and 0.870 against the validation set (Figure 5D). The AUC was even higher (0.918) against independent external validation data, the GSE151052 dataset (Figure 5E).

Figure 5 Identification of feature genes capable of diagnosing COPD. (A) Selection of optimal parameter (lambda) based on minimal criteria in the LASSO regression model. (B) LASSO coefficient profiles of 11 feature genes with non-zero coefficients. Receiver operating characteristic curves for feature genes applied to (C) the training data in the GSE47460 dataset, (D) the validation data in the GSE47460 dataset, or (E) the external validation data in the GSE151052 dataset. AUC, area under the receiver operating characteristic curve.

Immune Cell Infiltration in COPD

Since common genes were enriched in functions related to immune responses, we hypothesized that immune cell dysfunction might play a key role in COPD. Indeed, we found that infiltration by Th1 cells, follicular helper T cell (TFH), CD8+ T cells, and B cells was significantly greater in COPD patients than in controls (Figure 6A). In fact, these four immune cell types correlated positively with one another (Figure 6B). Feature genes correlated significantly with levels of immune cell infiltration in COPD patients (Figure 6C), with macrophages as the most abundant infiltrating cells (Figure 6D).

Figure 6 Differences in immune cell infiltration between COPD patients and controls. (A) Differences in immune cell infiltration between COPD patients and controls in the GSE151052, GSE38974, GSE47460, and GSE76925 datasets. Red indicates significant upregulation in patients; blue, significant downregulation. (B) Correlations and clusters among immune cell types in the GSE47460 dataset. (C) Correlations between immune cell infiltration and feature gene expression in COPD. *P < 0.05, **P < 0.01. (D) Levels of infiltration by 18 types of immune cells in COPD patients.

Regulation of Feature Genes by miRNAs

To identify how miRNAs may regulate feature genes, we started from 223 DEmiRs between COPD patients and controls in the GSE38974 dataset (Figure 7A). Then we predicted the downstream mRNAs targeted by the following six miRNAs with the largest log | (fold change) |: hsa-miR-1274a, hsa-miR-105-5p, hsa-miR-374a-3p, hsa-miR-422a-3p, hsa-miR-937-5p, and hsa-miR-923. Of the predicted mRNA targets, 162 were present among DEmRs modules C1_6, C1_3, and C1_2. These 162 downstream target genes were enriched in divalent inorganic cation homeostasis, cellular response to nitrogen compounds, as well as glutamate and glutamine metabolism (Figure 7B). Two of the genes, GPC4 and RS1, were predicted to be regulated by miR-374a-3p (Figure 7C). The AUC for miR-374a-3p was 0.89, suggesting the ability to diagnose COPD (Figure 7D). Finally, we generated a map of DEmiRs regulating DEmRs, which were involved in cell junction assembly, positive regulation of the MAPK cascade, and VEGFA−VEGFR2 signaling (Figure 7E).

Figure 7 Network of miRNAs that regulate COPD-associated mRNAs. (A) Differentially expressed miRNAs between COPD and controls in the GSE38974 dataset. Red dots indicate upregulated expression; green dots, downregulated expression. (B) The biological functions of target module genes enriched. (C) Predicted sites on the GPC4 and RS1 mRNAs where miR-374a-3p binds, based on Targetscan. (D) Receiver operating characteristic curve assessing the ability of miR-374a-3p to predict COPD. (E) Sankey map of miRNAs and the mRNAs that they regulate, together with the KEGG pathways in which they are involved.

Aberrant Gene Methylation in COPD

In the GSE63704 dataset, we identified 6495 DMPs between COPD and control samples, involving 4296 genes (Figure 8A). After retaining only DMPs that were increased or decreased between patients and controls in the opposite direction as common genes, we were left with 117 potential methylation markers (Figure 8B). Among them, we found that the feature gene PRRG4 was aberrantly methylated in COPD.

Figure 8 Aberrant methylation of common genes in COPD. (A) Differentially methylated probes (DMPs) between COPD patients and controls in the GSE63704 dataset, including hyper- and hypomethylated DMPs. (B) ∆Beta levels and expression levels of methylation markers.

Discussion

COPD is a common respiratory disease that seriously threatens human health and well-being, and it is an important public health problem.²⁶ Current treatments for COPD aim to improve symptoms and prevent exacerbations, but none is disease-modifying.²⁷ To provide the knowledge needed to improve therapy options, we compared gene expression profiles between COPD patients and healthy controls in order to identify coexpressed differentially expressed genes. MEGENA is an innovative method for analyzing coexpression networks that has advantages over weighted gene co-expression network analysis (WGCNA), facilitating efficient construction of large-scale coexpression plane filter networks while preserving gene interactions.²³ Our MEGENA identified 186 functional modules, the largest of which were C1_ 6, C1_ 3, and C1_ 2. Next we screened the differentially expressed genes using AUCs and LASSO regression to identify 11 feature genes, which we validated using an internal dataset (GSE47460) and external dataset (GSE151052). Feature genes have been widely used in multi-marker profiling studies.^28–30 Our feature genes in COPD showed high AUCs, suggesting their potential for clinical applications.

We found that many of the differentially expressed genes in COPD were associated with various immune responses. An increasing number of immune cell types have been associated with risk of COPD and prognosis of affected individuals.^31,32 In particular, COPD has been linked to elevated numbers of neutrophils, B cells, as well as CD4+ and CD8+ T lymphocytes in the lungs.³³ During adaptive immune responses in COPD patients, activated T cells can promote abnormal inflammatory responses and aggravate airway damage.³² Increased B cell counts in COPD patients have been associated with elevated IgA synthesis, which impairs mucosal immunity and may contribute to disease progression.³⁴ Consistent with this, our GSEA found that the intestinal immune network involved in IgA production was activated in COPD patients. Genetic alterations in leukocyte transendothelial migration pathways in smokers and COPD patients have been strongly associated with T cell levels and airway obstruction,^35,36 which our enrichment analyses support. In addition, asthma has been associated with the frequency and severity of COPD exacerbations,³⁷ which our analysis of DEmR enrichment supports.

Feature genes of COPD that we identified here have previously been linked to the disease, suggesting that our bioinformatics analysis is reliable. AQP1 is significantly up-regulated in COPD and has proven to be an effective therapeutic target.^38,39 Expression of CLEC4A is influenced by smoking,⁴⁰ which strongly correlates with risk of COPD and airflow obstruction.^41,42 GPC4, which remodels the extracellular matrix, is closely related to COPD disease regions.^43,44 GPR39 helps activate pro-inflammatory signaling pathways,⁴⁵ which may contribute to pathological inflammation in COPD. PDCL3 regulates expression of VEGF receptor 2, and it may promote pathological angiogenesis.^46,47 PGM1 has been linked to lung disease through multiple metabolic pathways.^48,49 PRRG4 has been shown to be differentially expressed in COPD patients,⁵⁰ and the present study not only confirms that finding but extends it by showing aberrant methylation of the gene. Therefore, the altered expression of PRRG4 in COPD may reflect altered gene methylation. ROR1, after binding to Wnt ligands, triggers non-canonical signaling cascades that increase the level of calcium or decrease the level of cGMP within the cell, and these changes are closely associated with COPD onset and progression.⁵¹ SEC14L4 is known to be differentially expressed in COPD.⁵² USP13, which inhibits autophagy, may be involved in COPD pathogenesis.^53,54

Among the 11 feature genes, RS1 has not previously been linked to COPD, and our study justifies further exploration of this potential link. In addition, our results suggest that the gene, together with GPC4, is regulated by miR-374a-3p. This miRNA has been shown to regulate inflammatory responses.^55,56 Further work should examine whether miR-374a-3p is involved in COPD via its ability to regulate the expression of RS1 and GPC4.

Our findings should be interpreted with caution in light of several limitations. We were limited to data available in public databases, and we did not validate key results using independent clinical samples or biochemical experiments. In addition, relevant clinical and follow-up data were unavailable for many of the samples that we analyzed, preventing us from examining potential relationships between feature genes and comorbidities or prognosis. Further study with larger samples and detailed follow-up should verify and extend the clinical utility of the feature genes that we identified here. Experimental studies in vivo and in vitro should explore how feature genes in COPD are regulated.

Conclusion

Many genes that are abnormally expressed in COPD are involved in immune responses, and we identified several feature genes that may be potential markers and therapeutic targets in COPD. In this way, our bioinformatics study generates numerous leads to guide future research into the disease and its treatment.