Molecular classification based on hypoxia-associated genes and construction of the prognostic model in Chronic Obstructive Pulmonary Disease

1 Introduction

Chronic obstructive pulmonary disease (COPD) is a pervasive and troublesome obstructive airway disease propelled by chronic inflammation of the respiratory tract that places an increasing burden on health care in industrialized and developing countries (Barnes, 2007; López-Campos et al., 2016). It is characterized by nonreversible airflow obstruction; primarily drags peripheral airways, that spells from air stagnation to dynamic hyperinflation and intermittent exacerbations (Tangedal et al., 2019; Fleming, 2019). COPD ails approximately 10 % of adults with age over 40 years old; and is a dominant count for admission, ranking as the third leading cause of death worldwide (Barnes et al., 2015). Although COPD typically presents as part of a multimorbidity pattern in old age; evidence has piled up to illustrate that some early life events disturb lung function, indicating that risk factors beyond those extensively reputed (examples abound: cigarette smoking and particles and gases inhaling) are glittering in the etiology of the disease (Svanes et al., 2010; Postma et al., 2015). Indeed, there is a long-standing interest in disentangling the genetic determinants from a wide range of COPD susceptibilities (Silverman, 2020). COPD is a heterogeneous disease that is subjected to genetic and environmental factors operating in a developmental climate (Decramer et al., 2012; Zhou et al., 2022). Painting a clear picture of genetic determinants in COPD allows for an unbiased evaluation of the crucial molecular determinants of disease pathobiology; which may yield encouraging insights into COPD pathogenesis (Huang et al., 2022).

COPD is frequently underdiagnosed or misdiagnosed. In general, spirometry is liable for the diagnosis of COPD; a post-bronchodilator FEV1/FVC < 0.70 acknowledges the existence of persistent airflow limitation and distinguishes patients with COPD suffering appropriate symptoms and predisposing risks (Buist et al., 2007; Singh, 2019). According to an investigation of the National Health and Nutrition Examination Survey (NHANES), more than 70 % of participants with chronic airway obstruction were not formally diagnosed with COPD (Martinez et al., 2015). The results of another analysis in five health plans have manifested that 32 % of patients with a new presumptive COPD diagnosis had been performed with spirometry to verify the detection (Han, 2007). Additionally, standardized dyspnea and symptom evaluation tools, for instance, the modified Medical Research Council (mMRC) dyspnea scale and COPD assessment test (CAT), could be capitalized on the stratification and surveillance of COPD progression (Mahler and Wells, 1988; Jones et al., 2009; Gupta et al., 2014). In the past 20 years, novel inhaled and oral medications (e.g., inhaled corticosteroids, ICS; and long-acting β2-agonists, LABAs) as well as emerging surgical and bronchoscopic procedures have been introduced for COPD treatment (Halpin, et al., 2020; Ferrera et al., 2021). Available pharmaceutical treatment for COPD can provide relief of burdensome symptoms and grind down exacerbation risk (Guo et al., 2023; Zhou et al., 2022; Lu, 2024).

In addition to medical therapy, supplementing oxygen is a momentous strategy of COPD interventions to consider in the context of severe hypoxemia. Severe hypoxia is a characteristic feature of COPD and is instrumental in inciting disease progression. Simultaneously, hypoxia is a ubiquitous attribute of a solid tumor microenvironment (TME) and is reputed as a pivotal factor in the hallmarks of cancers (Cui et al., 2022; Xu et al., 2020).

Recently, considerable research efforts are dedicated to unveiling the role of hypoxia in tumor development; however, only a few studies are in the running to dissect the function of hypoxia-related genes (HRGs) in COPD progression. A more comprehensive molecular classification and biological understanding of the HRGs in COPD could help to stratify patients into distinct phenotypes, thus contributing to the diagnosis and innovative drug development of COPD. In this study, we were posed to tap into the exact role and underlying mechanisms of HRGs in the pathological process of COPD. Notably, we interrogated the potential drug candidates targeting HRDEGs. The immune characteristics of COPD were also elucidated. Our research enlarged the current knowledge concerning the association between hypoxia and COPD.

2 Material and methods

3 Results

3.1

3.1 The acquisition of HRGs and DEGs in COPD

A total of 477 DEGs associated with COPD were discovered in the GSE54837 dataset (threshold as |log FC|＞ 1 and p-value < 0.05). The circular heatmap of DEGs between COPD and normal groups was visualized in Fig. 1A. 200 HRGs were identified by the hallmark hypoxias gene set. 6 overlapping genes were identified between DEGs and HRGs by Venn analysis (Fig. 1B). Subsequently, the gene expression profiles of 6 common HRGs were further exemplified in GSE54837. TPD52 and RORA were down-regulated in COPD samples while PPP1R15A, DUSP1, PDK3, and DDIT3 were up-regulated in COPD samples (Fig. 1C). Similarly, we discovered that PPP1R15A, DUSP1, and PDK3 were also up-regulated in the GES151052 dataset (Fig. 1D).

Close

Fig. 1

The acquisition of HRGs and DEGs in COPD. (A) The circular heatmap of DEGs in COPD groups in comparison with normal groups. (B) The Venn result of DEGs and HRGs. (C-D) The expression profile of 6 overlapping genes in GSE54837 and GSE151052.

Export to PPT

3.2

3.2 The identification of three imperative HRDEGs by machine learning method

We used RF, GLM, and SVM, three machine learning methods to figure out the HRGs that contributed greatly to COPD. The residual of the GLM method was the lowest of the three algorithms (Fig. 2A and 2B). The ROC analysis results showed that the AUC value of the GLM method was higher than RF and SAM (Fig. 2C). These findings demonstrated that the GLM algorithm performed excellently in the construction of the prognostic model of COPD. Results suggested that the importance of 6 overlapping HRDEGs varied according to three methods and they were different from each other (Fig. 2D). According to the GLM investigation result, the top three HRDEGs (TPD52, DDIT3, and DUSP1) were figured out to further analysis. ROC analysis indicated that TPD52 obtained the highest AUC value while DDIT3 obtained the lowest diagnostic value in COPD (Fig. 2E). The clinical nomogram of the three top HRDEGs was shown in Fig. 2F. Three HRDEGs could well distinguish COPD from normal samples. The calibration plot and ROC curves validated their good predictive power in COPD (Fig. 2G and 2H).

Close

Fig. 2

Identification of three crucial HRDEGs. (A) The reverse cumulative distribution of residual of three learning methods. (B) the boxes plot of residual for RF, GLM, and SVM. (C) The ROC curve of residual for RF, GLM, and SVM. (D) Top 6 HRDEGs lists by RF, GLM, and SVM screening. (E) The ROC curve of three HRDEGs identified by GLM. (F) The nomogram of three top HRDEGs. (G-H) The calibration plot and ROC analysis of prognosis nomogram.

Export to PPT

3.3

3.3 Potential drugs candidate of three HRDEGs

The screening of the TCMSP herbal database hinted at a natural compound targeting DUSP1, 17-beta-estradiol, which was the active ingredient of Amygdalus Communis Vas (also named bitter almond) (Fig. 3A). The protein product of DUSP1 could bind with 17-beta-estradiol by molecular docking analysis (Fig. 3B), indicating that Amygdalus Communis Vas could recognize DUSP1. We retained the top ten compounds of the combined score of the three HRGs in the DSigDB database (Fig. 3C). Molecular dock analysis found that diclofenac (top 1 compound) was able to simultaneously target TPD52 (Fig. 3D), DDIT3 (Fig. 3E), and DUSP1 (Fig. 3F). These findings revealed that diclofenac may be a candidate drug for treating hypoxia associated with COPD.

Close

Fig. 3

The drug-HRDEGs interaction analysis. (A) The interaction between 17-beta- estradiol and DUSP1. (B) Molecular dock results of DUSP1 with 17-beta-estradiol. (C) DSigDB database showed the top ten compounds of the combined score of three HRDEGs. (D-F) Molecular dock results of three HRDEGs and diclofenac.

Export to PPT

3.4

3.4 COPD subtypes analysis based on three HDRGs in GSE54837

The COPD samples were clustered by NMF, and it was found that the sample typing did the best performance at r = 3 (Fig. 4A). The matrix heat map also showed that COPD samples could be divided into three subtypes (named C1, C2, and C3) (Fig. 4B). The classification was verified by the NTP algorithm according to the top 100 upregulated genes of each subtype. The samples could also be broken down into three subtypes (Fig. 4C). Calibration plots displayed that NTP classification and NMF classification had a good correlation (r = 0.64, Fig. 4D). Among the three HDRGs, TPD52 was highly expressed in the C1 subtype (Fig. 4E), DDIT3 was predominantly elevated in the C3 subtype (Fig. 4F), and DUSP1 showed significant upregulation in the C2 subtype (Fig. 4G). Hypoxia score was higher in C2 and C3 subtypes than C1 subtype (Fig. 4H). Therefore, we defined C2 and C3 subtypes as the severe hypoxia subtype and C1 as the general hypoxia subtype. All three HDRGs had excellent predictive performance in predicting both hypoxia subtypes, with TPD52 being the most accurate (Fig. 4I).

Close

Fig. 4

Identification of COPD subtypes based on three HRDEGs (A) The sample typing result had a good performance of NMF at r = 3. (B) The heatmap of three COPD subtypes by NMF. (C) The heatmap of three subtypes was validated by the NTP algorithm. (D) The correlation analysis between NTP and NMF typing. (E-G) The expression of three HRDEGs in three COPD subtypes. (H) Hypoxia flux enrichment scores of three subtypes. (I) ROC curves of three HRDEGs in predicting COPD hypoxia subtypes.

Export to PPT

3.5

3.5 The endorsement of COPD typing based on three HRDEGs in GSE151052

We investigated the expression level of three HRDEGs in GSE151052 and COPD samples could be also classified into three subtypes (C1, C2, and C3) by the NMF typing method (Fig. 5A). Among the three HDRGs, TPD52 was highly expressed in C2 (Fig. 5B), DDIT3 was with highest expression in C1 (Fig. 5C), and DUSP1 expression was highest in C3 (Fig. 5D). The hypoxia score was higher in C1 and C3 subtypes (Fig. 5E). Similarly, we defined C1 and C3 as the severe hypoxia subtype and C2 as the general hypoxia subtype. TPD52 and DUSP1 predicted these two hypoxic subtypes well, and TPD52 was the most accurate prediction (Fig. 5F). This is essentially the same situation as in GSE54837.

Close

Fig. 5

Validation of COPD hypoxia typing in GSE15102 (A)The heatmap of three subtypes identified by three HRDEGs in GSE151052 via NMF. (B-D) The expression of three HRDEGs in three subtypes in GSE151052. (E) Hypoxia flux enrichment scores of three subtypes. (F) ROC curves of three HRDEGs in distinguishing COPD hypoxia subtypes.

Export to PPT

3.6

3.6 Immune characteristics analysis in COPD based on GSE54837

Among the three HRDEGs, TPD52 was positively correlated with immune features compared to DDIT3 and DUSP1, while DDIT3 and DUSP1 were negatively correlated with immune features (Fig. 6A). Relative to healthy individuals, immune features such as APC co-inhibition, type I IFN response, type II IFN response, and mast cells were more highly expressed in COPD (Fig. 6B). Immune cells such as B cells, CD8⁺ T cells, neutrophils, TIL, and immune signatures including APC co-stimulation, inflammation-promoting, MHC class I, type II IFN response, and other immune signatures differed in three subtypes (Fig. 6C). We found that most of the immune signatures with significant expression were highly expressed in the C1 subtype of COPD (Fig. 6D). It was concluded that patients with the general hypoxic subtype had better immunogenicity.

Close

Fig. 6

Immune cell and immune feature analysis of three HRDEGs in COPD dataset GSE54837. (A) The correlation between three HRDEGs and 16 immune cells and 13 immune signatures. (B) The expression profile of 16 immune cells and 13 immune signatures in COPD and normal groups. (C) The expression profile of 16 immune cells and 13 immune signatures in three subtypes. (D) Immune feature landscape in COPD subtypes.

Export to PPT

3.7

3.7 KEGG analysis in COPD hypoxia subtypes

We investigated the pathways enrichment between three COPD subtypes. Results indicated that the adrenergic angiotensin system and sulfur metabolism pathways were downregulated in C1 while the cytosolic DNA sensing pathway was downregulated in C2. Glycosaminoglycan biosynthesis keratan sulfate was down-regulated in C3 (Fig. 7A). Ribosome and protein export were up-regulated in C1 while porphyrin and chlorophyll metabolism pathways were up-regulated in C2 while non-homologous end- joining was up-regulated in C3 (Fig. 7B). This research uncovered that COPD with different hypoxia contexts exhibited metabolism pathways.

Close

Fig. 7

KEGG analysis in COPD hypoxia subtypes. (A-B) The KEGG enrichment results for COPD subtypes in the GSE54837 and GSE151052 datasets.

Export to PPT

3.8

3.8 Compounds sensitive to COPD with severe hypoxic context

The crucial metabolism signal, cAMP, was used to explore the association between COPD subtypes and small molecule compounds. TTNPB, arachidonyltrifluoromethane, MK.886, MS.275, and exisulind were found to be sensitive to COPD patients with severe hypoxic subtypes (Fig. 8). This could guide the treatment of patients with severe hypoxic subtypes of COPD.

Close

Fig. 8

Compounds sensitive to COPD with severe hypoxia.

Export to PPT

4 Discussion

Despite the yawning clout of cigarette smoking, genetic factors are also essential determinants of COPD. Indeed, numerous genomic regions that affect COPD susceptibility have been defined during the genome-wide correlation analysis, of which a quintessential example would be alpha-1 antitrypsin deficiency, a pervasive single-gene disorder among Northern Europeans and North Americans (Chiuchiolo and Crystal, 2016). Simultaneously, some genomic regions correlated with COPD-related phenotypes, for instance, quantitative measures of emphysema, have been identified (Silverman, 2018; Chen et al., 2019). Disentangling the functional variants and core genes within these linked regions is still riddled with challenges. The current study was designed to come to grips with the role of HRGs in COPD progression and their diagnostic value (He et al., 2019; Peng et al., 2022).

We first gleaned 477 DEGs between the healthy and COPD samples of the GSE54837 dataset. Further, 200 HRGs were obtained and 6 overlapping HRDEGs (TPD52, RORA, PPP1R15A, DDIT3, DUSP1, and PDK3) were screened between the DEGs and HRGs. Interestingly, the expression levels of TPD52 and RORA were tail off while counterparts of PPP1R15A, DDIT3, DUSP1, and PDK3 exhibited a mounting trajectory in COPD. The trend of PPP1R15A, DUSP1 and PDK3 were also observed in the validated set (GSE151052). PPP1R15A, serving as a stress-induced regulatory subunit, was sufficient to negatively modulate eIF2α dephosphorylation (Hölzer et al., 2016). Phosphorylation of eIF2α is an evolutionarily conserved and critical cellular defense system that can counteract multiple forms of stress (Carrara et al., 2017). Studies have shown that PPP1R15A was sinking in lung fibroblasts of patients with idiopathic pulmonary fibrosis (IPF; an unrelenting; chronic, and progressive lung disease that could be triggered by smoking) (Monkley et al., 2021; Richeldi et al., 2017). TGF-β-educated fibroblasts further shrunk PPPIR15A. The absence of PPPIR15A or its pharmacological handcuff aggravated experimental lung fibrosis. Whether such an effect is thrown on COPD remains currently unclear.

Three machine learning algorithms (GLM, RF, and SVM) demonstrated the varied clout of 6 HRDEGs in COPD progression; and, consequent to that, 3 paramount HRDEGs (TPD52, DDIT3, and DUSP1) were finally identified based on GLM analysis as it presented the largest area under ROC curve. Taken as a whole, the top three key HRDEGs, albeit with distinct predictive probability, performed well in distinguishing patients with COPD. The tumor protein D52 (TPD52) family was generally interpreted as a crucial driver in the proliferation and metastasis of a wide range of cancers (Fan et al., 2021; Byrne et al., 2005). TPD52L2, a member belonging to the TPD52-like protein family, was overexpressed and associated with unfavorable prognosis in a variety of cancers, including lung adenocarcinoma, glioblastoma, and prostate cancer (Zhong et al., 2021; Qiang et al., 2018; Ren et al., 2017). However, studies centered around TPD52 in COPD are still scarce. Our study gives a window to decode the role of TPD52 in COPD progression and its potential in clinical diagnosis.

In recent years, a large number of studies concerning traditional herbal medicine have been published, and many attractive applications have been explored (Li and Weng, 2017; Zhao et al., 2019). Having a series of surfing among the TCMSP herbal medicine, it was observed that, 17-β- estradiol, the biologically active component of Amygdalus Communis Vas, could bind to DUSP1. Low-expressed DUSP1 galvanized alternation later expression pattern of genes engaged in certain biological pathways, including angiogenesis, cellular communication, and tyrosine-kinase receptor activity (Moncho-Amor et al., 2011). Such alternation was largely boiled down to the administration of c-Jun-N-terminal kinase and/or p38 activity by DUSP1 (Hao et al., 2015). It has been reported that upregulation of DUSP1 abundance and consequently altered patterns of JNK and p38 activities are strongly associated with NSCLC initiation, invasion, and metastasis (Moncho-Amor et al., 2011). A proportion of patients with COPD progressed rapidly to SCLC, despite the molecular mechanisms linking COPD with lung cancer development being far from clear (Cho et al., 2018; Xiao et al., 2017). Overall, we identified a novel kind of herbal medicine that can target DUSP1, thus manipulating its expression levels. Moreover, we found that diclofenac, with the highest combined score, could act on TPD52, DDIT3, and DUSP1.

Given the importance of HRGs in the progression of COPD, it could be advantageous to develop key HRGs as biomarkers that contribute to distinguishing COPD patients. We defined three hypoxia subtypes (namely C1, C2, and C3) by consensus clustering based on the expression landscape of three HRDEGs. In GSE54837, TPD52, DDIT3, and DUSP1 showed an overexpression tendency in the C1, C3, and C2 subtypes, respectively. According to the hypoxia score, we nominated C2 and C3 as the severe hypoxia subgroup and C1 as the normal hypoxia subgroup. The three pivotal HRDEGs, in particular, TPD52, performed well in stratifying COPD patients into the two subgroups with varying degrees of hypoxia. Such characteristics of the HRDEGs exhibited a high diagnostic value in terms of COPD. In the validated set GSE151052, we found that the expression of TPD52 was upregulated in C2, and DDIT3 was overexpressed in C3. We defined C1 and C3 as the severe hypoxia subgroup and C2 as the normal hypoxia subgroup. Within this framework, TPD52 also showed considerable precision in identifying the specific hypoxia subgroups among COPD patients, reinforcing its utility in the molecular stratification of the disease based on hypoxia status.

There are, to date, very limited studies comparing the pathology between COPD and other causes of disease other than cigarette smoking (e.g., immunopathology), to this end, we were in the running to tackle this issue through ascertaining immune landscape among different subtypes (Birring et al., 2002; Rivera et al., 2008). Interestingly, TPD52 exhibited a positive correlation with immune characteristics, while DDIT3 and DUSP1 displayed a negative association. Immune features like APC co-inhibition, type-I IFN response, type II IFN response, and mast cells contributed to the identification just mentioned. These could be employed to estimate the immunity of patients with COPD. Alternatively, hypoxia subtypes could be cooperated with the immune features to make a comprehensive consideration of the patient’s condition. Such a combination might not only identify patients who could benefit from immunotherapy but also cast light on predicting the prognosis of COPD patients.

Enrichment analysis demonstrated the difference in pathways among the three subtypes, further confirming the specific immune landscape in them. Aside from these, we also explored some agents, for example, TTNPB, MK.886, and MS.275, that were engineered to alleviate the progression of COPD patients in the severe subtype.

5 Conclusion

In conclusion, our study emphasized the role of HRGs in COPD progression and identified three COPD hypoxia subtypes based on predominant HRDEGs (TPD52, DUSP1, and DDIT3). The immune characteristics of COPD hypoxia subtypes differed from each other. These observations may provide clues for the diagnosis and treatment of patients with COPD.

Availability of data and materials: The dataset supporting the conclusions of this article is included within the article.

Funding

None.

CRediT authorship contribution statement

Zhongshuai Fu, Dongsheng Song conceived the study and wrote the article. Qingrong Cui, Danbo Li, Beilei Wang contributed to data interpretation. Qingwei Zhou and Xianfei Ding revised the article.