Critical Biological Functions and Clinical Implications of Epigenetic-Related Candidate Biomarkers in Chronic Obstructive Pulmonary Disease: Integrated Machine Learning Screening and Basic Experimental Validation

Introduction

Chronic obstructive pulmonary disease (COPD) is a complex condition marked by enduring and advancing restriction of airflow. It has emerged as the third most prevalent cause of mortality globally and is projected to lead to more than 5 million fatalities by 2060.¹ With the development of COPD management, the symptoms and further risks of COPD patients have been significantly improved.² However, few preventative and treatment are available to reduce disease progression and improve survival period and life quality. The existing biomarkers for COPD often lack sufficient sensitivity for early diagnosis and are unable to reflect the complex pathological mechanism of the disease.³ Hence, it is an urgent need for further exploration of novel therapeutic targets in COPD and development of biomarkers with clinical utility.

Epigenetic modifications, such as DNA and RNA methylation, play a critical role in regulating gene transcription programs and ultimately determining cell fate. Recent research on epigenetic modifications has not only shed light on the various regulatory roles of nucleic acids within living organisms but has also showcased their promising uses in clinical practice. Advanced methods, such as loci-specific, genome-wide methylation analysis and single-cell epigenomics, have been employed to detect epigenetic abnormalities in various contexts, including early diagnosis of malignancies, disease prognosis, evaluation of drug responses, and monitoring of complications.^4–6 Several epigenome-wide association studies (EWASs) have identified cytosine-phosphate-guanine (CpG) probes related to pulmonary function.^7,8 However, most EWASs faced limitations, such as reliance on blood samples and being conducted predominantly in populations of European ancestry. Furthermore, the underlying mechanisms of epigenetic modifications in COPD remain poorly understood due to the absence of a stable short-term model. Exposure to cigarette smoke (CS) in mouse lungs and human lung epithelial cells has been demonstrated to induce a pro-inflammatory gene expression pattern via histone acetylation.⁹ Emerging evidence suggests that noncoding RNAs, including miRNAs and lncRNAs, may be key regulators of epigenetic and gene transcription. When exposed to cigarette smoke extract, human bronchial epithelial cells exhibited upregulation of miR-101 and miR-144.¹⁰ However, CS-induced models may not fully reflect the in vivo functionality of epigenetic modifications in COPD, as COPD is a heterogeneous lung condition arising from gene-environment interactions.

Epigenetic-related genes (ERGs) were defined as a set of genes encoding core epigenetic regulators (DNA/histone mark readers/writers/erasers/remodellers).¹¹ To ensure the reliability and comprehensiveness of the ERG set, we retrieved candidate genes from the EpiFactors database:¹² a well-recognized, manually curated database of human epigenetic factors and their complexes, which has been extensively applied to dissect the biological functions of epigenetic modifications in various pathological processes. Various types of machine learning algorithms, such as supervised, unsupervised, semi-supervised and reinforcement learning, in addition to deep learning, which is a member of the machine learning family, can intelligently and comprehensively analyze large amounts of data.¹³ In this study, we used machine learning algorithms, including the least absolute shrinkage and selection operator (LASSO), Support Vector Machine Recursive Feature Elimination (SVM-RFE), and Boruta analysis to screen the relatively stable and important characteristic genes in COPD. To explore the underlying the regulatory mechanisms, we performed an immune infiltration analysis and constructed potential regulatory network based on the biomarkers. Finally, we investigated the clinical relevance of biomarkers, assessing their utility in disease prediction and as potential drug targets. In general, this study has filled a crucial gap in the field of epigenetics research on COPD. By conducting a comprehensive network-level analysis of ERGs, it addressed the absence of such systematic research and overcame the limitations of traditional statistical methods through machine learning, thereby identifying reliable biomarkers. Moreover, this study bridged the gap between the discovery of epigenetic biomarkers and their functional interpretation, laying a theoretical foundation for translating these insights into clinical diagnosis and treatment of COPD.

Materials and Methods

Data Acquisition

A total of three COPD transcriptome microarray datasets (GSE19407, GSE11784 and GSE20257) were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). All three datasets were composed of the transcriptome of small airway epithelium in healthy non-smokers, healthy smokers and smokers with COPD, with the GSE19407 dataset consisting of 105 normal and 22 COPD samples, the GSE11784 dataset composed of 135 normal and 22 COPD samples, and the GSE20257 dataset consisting of 112 normal and 23 COPD samples. The 720 ERGs were downloaded from the EpiFactors database (http://EpiFactors.autosome.ru).

Identification of Differentially Expressed Epigenetic-Related Genes (DE-ERGs) in COPD

Differentially expressed genes (DEGs1, DEGs2, and DEGs3) in the COPD group and normal group were identified in three datasets (GSE19407, GSE11784, and GSE20257) utilizing the “limma” package (version 3.54.0). The statistical significance criterion was |log₂FC| > 0.5, and P < 0.05. Volcano maps of all in 3 datasets were generated utilizing the “ggplot2” package (version 3.4.1), and the “pheatmap” package to generate heat maps. In order to obtain the common DEGs, the up-regulated genes and down-regulated genes were intersected, respectively, in the three datasets. Then DEGs and 720 ERGs were further intersected to obtain the DE-ERGs in COPD. The “ggplot2” software package was used to create Venn diagrams.

Enrichment Analysis and Protein-Protein Interaction (PPI) Network Analysis

Gene ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) were utilized to analyze the role of DE-ERGs in biological processes and the biological pathways. GO and KEGG enrichment analyses of DE-ERGs were conducted utilizing the R package “clusterProfiler” (version 4.2.2) and to look for common functions and related pathways between genes, respectively. The main enriched functions and pathways of DE-ERGs were screened at P < 0.05, and the GO-enriched entries and KEGG-enriched entries were visualized using the R language “enrichplot” package (version 1.18.3) and “ggnewscale” package, respectively. The PPI network of DE-ERGs was carried out utilizing the STRING database (http://string.em bl.de/) and visualized utilizing Cytoscape software, which utilized the plug-in cytoHubba to detect node genes in the PPI network (medium confidence = 0.400).

Screening for Biomarkers Using Three Machine Learning Algorithms

Based on the nodal genes screened by the PPI network, we used the “glmnet”, “e1071” and “Boruta” packages for the LASSO, SVM-RFE, and Boruta analysis, respectively. The LASSO model adopted 10-fold cross-validation with α fixed at 1 and finally selected the model corresponding to the maximum λ value (0.0067) within one standard deviation of the minimum cross-validation deviance. SVM-RFE employed the default radial basis kernel function. During the iterative feature elimination process, the importance of features was evaluated through internal cross-validation. After 5-fold external cross-validation and repeating this process 5 times to obtain a stable generalization error estimate, the feature genes corresponding to the minimum average error are selected as the final result. Boruta employed default parameters (ntree = 500, maxRuns = 100), using the expression levels of candidate genes as the feature and the sample grouping as the response variable for feature selection. It screened out the important feature genes by comparing the importance Z scores of the original features with those of the shadow features. Subsequently, the common genes identified by the three machine learning methods were selected as a biomarker. Based on the datasets GSE19407, GSE11784, and GSE20257, ROC curve analyses were carried out on the biomarkers utilizing the “pROC” software package, and the prediction efficiency of biomarkers was verified according to the area under the curve (AUC) values. The R package “ggplot2” was further used to plot the expression box plots of biomarkers in the three datasets.

Chromosomal Localization Analysis and GeneMANIA Analysis of Biomarkers

Human chromosome data (UCSC.HG38.Human.CytoBandIdeogram) were imported using the R package “Circos”, and circle maps were drawn to visualize the distribution of biomarkers on chromosomes. Using the online tools on the GeneMANIA website (https://genemania.org/), we predicted the genes associated with biomarkers function, and constructed gene-gene co-expression networks for the biomarkers.

Gene Set Enrichment Analysis (GSEA) Enrichment and Immunoinfiltration Analysis

The genome named “c2.cp.kegg.v7.4.symbols.gmt” was retrieved from the molecular characterization database, and GSEA was performed on the biomarkers to elucidate the significant functional and difference pathway of the biomarkers, using the significance criterion of P < 0.05. The “ssGSEA” package was utilized to evaluate the count of immune cell infiltrations in the GSE19407 dataset, which was then visualized through a stacked plot. Subsequently, an analysis was conducted to compare the differences in each immune cell type between the COPD and control groups (adj.P < 0.05), and expression box plots were drawn using the “ggplot2” package. Finally, Spearman analysis was utilized to observe the association between biomarkers and differentially immune cells.

Construction of Regulatory Network for COPD

We used the staBase (https://rnasysu.com/encori/) database to predict microRNAs (miRNAs) targeting biomarkers, followed by the ChEA3 database (https://maayanlab.cloud/chea3/) to predict transcription factors (TFs) that could regulate biomarkers, and Cytoscape was utilized to create a TF-miRNA-mRNA regulatory network based on the predicted TFs, miRNAs and mRNAs.

Construction of Gene-Disease Network and Gene-Drug Prediction

The DisGeNET database (https://www.disgenet.org/) collected and integrated disease-gene association information from multiple data sources, through which biomarkers-disease association networks were derived. Based on the biomarkers of COPD, the CTD database (https://ctdbase.org/) was utilized to predict potential drug components for the treatment of COPD. The Cytoscape software (version 3.8.2) was utilized for visualizing the gene-drug networks.

Expression of Biomarkers in Clinical Samples

Source of Human Samples

Human induced sputum samples were gathered from a cohort of COPD patients as well as from healthy human volunteers without known cardiopulmonary disease, following the protocol approved by the institutional review board at Zhujiang Hospital (Approval Number: 2023-KY-160-02). Consistent with the global initiative for Chronic Obstructive Lung Disease (GOLD) 2023, the study enrolled participants included individuals aged 40–60 years who had a specialist diagnosis of stable COPD. All participants or their legally authorized representatives provided written informed consent. The control group comprised individuals aged between 40 and 55 years who had no prior history of COPD and had a baseline FEV1/FVC ratio of 0.7 or higher. The characterization details of these participants were shown in Table 1. Induced sputum was conducted with nebulized hypertonic 3% saline for 20 minutes. Albuterol (180 mg via metered dose inhaler) was administered 15 minutes before sputum induction. The sputum samples were dispersed in an equal volume of 0.1% DTT and the supernatant removed after the sample was centrifuged at 1000 × g. The cell pellet was stored at −80°C for subsequent analysis.

Table 1 Characterization Details of Participants

Quantitative Reverse Transcriptase Polymerase Chain Reaction (qRT-PCR) Analysis

Trizol (Invitrogen, USA) was utilized to extract total RNA, cDNA was synthesized from 1 µg of RNA extracted from sputum samples using the Prime Script RT Kit (Takara, Dalian, China). Real-time PCR (RT-PCR) was carried out utilizing Bio-Rad CFX Connect detector system with SYBR Green universal PCR mix (Takara). Primers were designed and acquired from Sangon Biotech (Shanghai, China). UBE2T, HMGN4, TLE1, CIT and TFPT mRNA expression was normalized to 18S expression.

Fold changes in mRNA expression were calculated utilizing the 2^−ΔΔCt method, The following primers were displayed in Table 2.

Table 2 The Primer Sequences of Biomarkers

Statistical Analysis

All analyses were performed utilizing R version 4.2.1. The Wilcoxon test was utilized for significance analysis to compare the differences between the two groups of samples. The data were presented as mean ± standard deviation. For comparisons between two independent groups, normally distributed continuous data were analyzed using an unpaired Student’s t-test, whereas non-normally distributed data were assessed with the Mann–Whitney U-test. Normality was verified using the Shapiro–Wilk test (or visual inspection of Q–Q plots, as appropriate). All statistical analyses were performed using GraphPad Prism software. The statistical significance was demonstrated at a level of P < 0.05.

Results

Determination of DE-ERGs in COPD

The GSE19407, GSE11784, and GSE20257 datasets yielded 2466 DEGs1 (including 2250 up-regulated and 216 down-regulated genes), 2071 DEGs2 (including 1829 up-regulated and 242 down-regulated genes), and 2347 DEGs3 (including 2143 up-regulated and 204 down-regulated genes), respectively. The volcano maps and heat maps showed the differential expression of genes in the three datasets, respectively (Figure 1A and B). A total of 1677 DEGs in COPD were obtained by taking the intersection and then the concatenation of the up-regulated and down-regulated genes of DEGs1, DEGs2, and DEGs3, respectively (Figure 1C). Then, 1677 DEGs and 720 ERGs were taken to intersect to obtain 23 DE-ERGs in COPD (Figure 1D).

Figure 1 Identification of differentially expressed genes (DEGs) in the GSE19407, GSE11784, and GSE20257 datasets. (A) Volcano plot of DEGs in 3 datasets; (B) Heatmap of DEGs in 3 datasets; (C) Venn diagram of up- and down-regulated DEGs in 3 datasets. (D) The overlapping genes between DEGs and ERGs.

Functional Enrichment and PPI Analysis of DE-ERGs

GO analysis was mainly enriched to the entries of regulation of single stranded viral RNA replication via double stranded DNA intermediate, deoxycytidine deaminase activity and viral RNA genome replication. The top 10 entries were visualized according to the p-value value from smallest to largest, and the results were shown in Figure 2A. KEGG analysis mainly enriched 4 pathways: HIV-1 viral life cycle, FOXO signaling pathway, viral carcinogenesis, and human immunodeficiency virus 1 infection (Figure 2B). The PPI network showed 13 node genes and 10 interacting pairs, with the strongest interactions for APOBEC3A and APOBEC3B (Figure 2C and Table S1).

Figure 2 Enrichment analysis of differentially expressed epigenetic related genes (DE-ERGs) and construction of protein-protein interaction (PPI) network. (A) Gene Ontology (GO) enrichment terms; (B) Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment pathways; (C) PPI network construction.

A Total of 5 Biomarkers of COPD Were Screened by Machine Learning Algorithms

Based on the GSE19407 dataset, the 13 node genes were analyzed by machine learning, and 13, 7 and 6 characteristic genes were obtained by SVM-RFE analysis, LASSO analysis, and Boruta analysis, respectively (Figure 3A–D). Subsequently, the genes derived from the three machine learning algorithms were taken as intersections to obtain 5 biomarkers (HMGN4, CIT, TLE1, TFPT, and UBE2T) (Figure 3E). The ROC curves indicated that the AUC values of the biomarkers in three datasets exceeded 0.7 (Figure 4A–C), suggesting that 5 biomarkers had acceptable diagnostic capabilities of COPD. The expression levels of 5 biomarkers were up-regulated in all three datasets (Figure 4D–F).

Figure 3 Construction of Machine Learning Diagnostic Models. (A-D) Feature genes screening in the (A) support vector machine recursive feature elimination (SVM-RFE); (B-C) least absolute shrinkage and selection operator (LASSO); and (D) Boruta algorithms. (E) Three machine learning algorithms selected feature genes and the Venn diagram of their intersection.

Figure 4 Receiver operating characteristic (ROC) analysis and expression validation of five key genes in the GSE19407, GSE11784, and GSE20257 datasets. (A-C) ROC curve of key genes in three datasets; (D-F) Validation of expression of key genes in three datasets. *** P < 0.001, **** P < 0.0001.

Chromosomal Localization and GeneMANIA Network Analysis of Biomarkers

The circle diagram showed the localization of the 5 biomarkers in the chromosomes, UBE2T, HMGN4, TLE1, CIT, and TFPT were localized on 1, 6, 9, 12, and 19 chromosomes, respectively (Figure 5A). The GeneMANIA network map showed a total of 20 genes associated with biomarker function, with HES1, FANCL and SRAI had the stronger interactions with these biomarkers (Figure 5B).

Figure 5 Chromosomal localization, GeneMANIA network, and Gene set enrichment analysis (GSEA) enrichment analysis of key genes. (A) Chromosomal distribution of key genes; (B) Construction of GeneMANIA Network for key genes; (C-G) GSEA analysis enrichment of key genes.

GSEA Analysis of 5 Biomarkers

After GSEA analysis of the 5 biomarkers, we found that the HMGN4 was enriched in total 84 pathways, mainly enriched in oxidative_phosphorylation, pentose and glucuronate interconversions, etc. CIT was enriched 58 pathways, mainly enriched in MAPK signaling pathway, ascorbate and aldarate_metabolism. TLE1 was enriched in 72 pathways, mainly in MAPK signaling pathway, ascorbate and aldarate metabolism, etc. TLE1 was enriched in 72 pathways, mainly in metabolism of xenobiotics by cytochrome p450, ascorbate and aldarate metabolism, etc. TFPT was enriched in 91 pathways, mainly in Huntingtons disease, oxidative_phosphorylation. TFPT was enriched in 91 pathways, mainly in Huntingtons disease, oxidative phosphorylation and other pathways. UBE2T was enriched in 81 pathways, mainly in ribosome, proteasome and other pathways. We selected the top 5 pathways for display (Figure 5C–G).

Immune Infiltration Analysis of Biomarkers

The stacked plot indicated that the count of 28 immune cells between the COPD and control groups in the GSE19407 dataset (Figure 6A). The box plot showed that the 11 immune infiltrating cells such as activated B cell, central memory CD8 T cell, effector memory CD4 T cell, etc. were significantly different between the groups (Figure 6B). Overall, the 5 candidate biomarkers (UBE2T, TFPT, TLE1, CIT, and HMGN4) exhibited predominantly positive correlations with most immune cells. Notably, HMGN4 showed no significant correlation with activated B cells, but had significant positive correlations with the other 10 types of immune cells. Among them, the correlations with type 2 T helper cells and effector memory CD4 T cells were the strongest (cor > 0.5, P < 0.05). Additionally, activated B cell had significant negative correlations with the other 4 candidate markers except for HMGN4 (Figure 6C). These results suggest that these 5 candidate biomarkers may participate in the immune disorder process of COPD by regulating the infiltration and activation status of immune cells and have potential immunoregulatory significance.

Figure 6 Immunoinfiltration analysis. (A) Percentage of immune cell infiltration in normal and COPD samples; (B) Differences in immune cell infiltration between normal and COPD groups; ns, not significant; * P < 0.05. (C) Correlation analysis of key genes and immune cells. Red indicates positive correlation, blue indicates negative correlation; results marked with “×” are not statistically significant.

Figure 6 continued.

TF-miRNA-mRNA Regulatory Network of Biomarkers

To further explore the regulatory mechanisms of biomarkers, miRNAs and TFs targeting biomarkers were predicted. We predicted 20, 9, 2, and 3 miRNA for CIT, HMGN4, TLE1, and UBE2T, respectively, but no miRNAs were predicted for the TFPT. The TF-miRNA-mRNA regulatory network showed 127 nodes with 167 action pairs, and CTCF could simultaneously regulate the expression of UBE2T and TLE1 (Figure 7A and Table S2).

Figure 7 Construction of Key Gene Related Networks. (A) Transcription factor (TF)-microRNA (miRNA)-mRNA regulatory network; (B) Gene disease association network; (C) Drug-gene interaction.

Figure 7 continued.

Gene-Disease Association Networks and Gene-Drug Prediction Networks for Biomarkers

The gene-disease network had 28 nodes and 25 edges, and we selected the TOP5 to construct the network and display it, and showed that both the TLE1 and the HMGN4 predicted neoplasms, and both HMGN4 and UBE2T predicted liver carcinoma (Figure 7B). The gene-drug network demonstrated 30 nodes and 37 action pairs, with 18, 6, 13, and 2 potential drug components predicted for the CIT, TFPT, TLE1, and UBE2T, respectively. Among the drugs associated with the development of COPD, bisphenol a and valproic acid were predicted by CIT, TFPT and TLE1, tobacco smoke pollution by both TLE1 and UBE2T, and furan by CIT and TFPT (Figure 7C).

Validation of the Expression for Biomarkers in Clinical Samples

To further verify the difference in the expression of biomarkers between normal and COPD samples, we collected clinical samples for qRT-PCR experiments. The results showed that the expression of 5 biomarkers (HMGN4, CIT, TLE1, TFPT, and UBE2T) were significantly higher in COPD samples than in normal samples, which were consistent with the analysis results in public databases, thus further improving the stability and reliability of our findings (Figure 8).

Figure 8 Validation of the expression for biomarkers in clinical samples. * P < 0.05; ** P < 0.01.

Discussion

With a global prevalence exceeding 10%, COPD ranks as the third most significant cause of mortality worldwide.¹⁴ Due to the lack of typical early symptoms and reliable biomarkers, once COPD is diagnosed, most patients are at intermediate or advanced stages, often missing the optimal treatment window. Recent study suggested that epigenetic modifications were associated with reduced lung function and inflammatory response.^15,16 Therefore, seeking for effective and novel biomarkers linked to epigenetic modification has considerable promise, potentially providing targets of diagnosis and therapy in COPD, thereby enhancing clinical outcomes.

In this study, common DEGs were found in three datasets of COPD patients and normal samples, with 1535 up-regulated genes and 142 down-regulated genes, respectively. Among these DEGs, 23 genes were associated with epigenetic modifications. To understand the biological function of 23 DE-ERGs, GO and KEGG analyses were performed. Using GO enrichment analysis, this study identified that DE-ERGs were primarily associated with “regulation of single stranded viral RNA replication via double stranded DNA intermediate”. This may suggest that COPD could increase the risk and worsen clinical outcomes of RNA virus infection. Cigarette smoke (CS) exposure, the leading cause of COPD, impairs antiviral responses and increased viral replication in bronchial epithelial cells.¹⁷ Exacerbations of the disease, marked by sudden decline in lung function, are related to RNA virus infections, such as rhinovirus and respiratory syncytial virus (RSV),¹⁸ contributing to COPD progression.¹⁹ The GO analysis also predicted the highest scoring molecular function as “deoxycytidine deaminase activity”.²⁰ Deoxycytidine kinase was found upregulated in lungs of COPD patients, possibly due to hypoxia, and may contribute to increased apoptosis.²¹ Detecting genes modifying deoxycytidine in sputum might serve as predictive and pharmacodynamic biomarkers for COPD patients.²² The present study further performed KEGG pathway enrichment analysis, which identified that the DE-ERGs were associated with “HIV-1 viral life cycle” and “FOXO signaling pathway”. HIV-1 infection serves as a significant risk factor of COPD. HIV-1 virus infection enhances the retention of CD8+ T cells within the airway mucosa and induces the expression of transforming growth factor-β1 (TGFβ1) in airway epithelial cells,^23,24 which might be implicated as a mechanism of COPD development. Alterations in the microbiota of various body sites in HIV patients with COPD may indirectly promote the progression of chronic bronchitis and emphysema.^25,26 Additionally, the FOXO signaling plays a critical role in extensive airway remodeling, a hallmark of COPD.²⁷ Activation of the FOXO signaling pathway of the proteasome system is involved in the atrophy of COPD peripheral muscle cells, primarily through enhanced autophagy and oxidative stress.^28,29 Inhibition of the FOXO signaling pathway has been shown to protect lung function decline and exert anti-pulmonary fibrosis effect.³⁰ Collectively, these enrichments suggest that epigenetic-related biomarkers may influence COPD progression through these pathways, though direct causal links remain to be validated.

Biomarkers associated with epigenetic modification in COPD were screened through bioinformatics analyses, including HMGN4, CIT, TLE1, TFPT and UBE2T. HMGN4 performs its biological functions by facilitating gene transcription through its binding with nucleosomes. Despite being discovered as a novel epigenetic regulator in 2001, its biological roles have been scarcely investigated.³¹ A recent study has revealed a positive association between HMGN4 and STAT3, which promoted cell proliferation and tumor progression.³² However, there are few studies of HMGN4 in COPD. In our research, we observed elevated expression of HMGN4 in COPD samples compared to normal controls, which suggests a potential role for HMGN4 in COPD progression. CIT, a crucial component of the midbody, is indispensable for cytokinesis.³³ As a Rho effector, CIT was shown to directly phosphorylate eNOS-Thr497, associated with vascular dysfunction.³⁴ Overexpression of CIT may potentially contribute to the increased risk of cardiovascular diseases, which are prevalent comorbidities in COPD patients. TLE1 acted as a corepressor by interacting with various DNA-binding transcription factors, inhibiting the nuclear factor-κB pathway linking with chronic inflammation,^35,36 and regulating the Wnt pathway inducing lung regeneration in COPD.^37,38 TFPT had been recognized as a molecular fusion partner of TCF3 and is implicated in processes like cell proliferation and the induction of programmed cell death.³⁹ Moreover, studies had shown that TCF3 was highly expressed in lung cancer, suggesting that TFPT played a role in lung tissue diseases.⁴⁰ UBE2T facilitated the repair of DNA interstrand crosslinks and stalled replication forks.⁴¹ Studies have identified UBE2T as a major oncogene in lung cancer,⁴² and its high expression in COPD samples might promote disease progression. In summary, the upregulation of these biomarkers in COPD likely reflects their potential involvement in driving disease processes, albeit further mechanistic studies are required.

To further explore the biological processes involved in biomarkers, GSEA analysis was performed, we observed increased enrichment of oxidative phosphorylation for HMGN4, TFPT and UBE2T in COPD patients, which was supported by previous studies demonstrating upregulated expression of the mitochondrial electron transport chain and oxidative phosphorylation in CS-induced mice and COPD patients.^43,44 Enhanced oxidative phosphorylation leads to increased ROS, important to aged-related comorbidities, such as lung cancer, bronchiectasis, periodontitis, diabetes, and so on.^45–47 These findings aligned with our results, suggesting that our discoveries were accurate as well as providing further evidence and a theoretical basis for COPD pathogenesis.

Studies have shown that COPD patients exhibited abnormal immune responses, which could lead to airway inflammation and lung damage in COPD.⁴⁸ So, immune infiltration analysis was performed, we identified 11 types immune infiltration cells in COPD samples and normal samples, including activated B cell, central memory CD8 T cell, effector memory CD4 T cell, and so on. Further analysis revealed that TLE1 and HMGN4 exhibited stronger positive correlation with effector memory CD4 T cell and immature B cell, and numerous immunohistochemical studies have shown increased volume fractions of CD4⁺ T cells, B cells in the small airways and parenchyma of COPD patients.^49,50 Furthermore, type 2 T helper cells showed a significant and strong positive correlation with all 5 candidate biomarkers, which was highly consistent with the conclusion from previous studies that the imbalance of T helper cell cytokines (including Th1/Th2 and Th17/T regulatory cells) is involved in the pathophysiological process of COPD.⁵¹ Additionally, HMGN4 showed a significant positive correlation with 10 immune cell subtypes other than activated B cells, including effector memory CD4 T cells and natural killer T cells. Yang et al’s research in hepatocellular carcinoma also reported a similar positive correlation between HMGN4 and T cells, further supporting the results of this study. These observational correlations suggest that the identified epigenetic-related candidate genes may be indirectly associated with the dysregulation of immune cell subpopulations in COPD, which warrants further mechanistic validation.

However, most current studies on epigenetics in COPD remain descriptive or retrospective, lacking systematic screening of epigenetic-related biomarkers using machine learning-based bioinformatic strategies with large-scale datasets. To address this gap, this study leveraged the GEO public database to identify candidate epigenetic biomarkers associated with lung function decline and inflammatory response in COPD through an integrated machine learning framework. These findings provide a foundation for future experimental validation and may inform the development of diagnostic or prognostic tools for COPD. However, several limitations should be acknowledged in the present study. However, this study is an exploratory research primarily based on bioinformatics analysis and preliminary experimental validation. The conclusions drawn only reflect correlations rather than definite causal relationships, and the underlying regulatory mechanisms as well as clinical translational value still require further clarification through more in-depth mechanistic investigations and large-sample validation in subsequent studies. Although the qRT-PCR verification in this study was based on a limited sample size, the results were statistically significant and consistent with those from public databases. Future studies will need to conduct further verification of gene expression levels and clinical applicability through larger-scale multi-center cohorts. In addition, the dataset from GEO database used in this study lacks information on pulmonary function, laboratory findings, and radiologic data, limiting the scope of our conclusions.

Conclusion

To sum up, this study screened five epigenetic-related candidate genes (HMGN4, CIT, TLE1, TFPT and UBE2T) for COPD using bioinformatics analysis, which were subsequently validated in a small cohort of COPD patients. Based on GEO database expression profiles and functional annotation, these genes were found to be associated with virus infections, DNA modifications, FOXO signaling, and oxidative phosphorylation. These candidate biomarkers might contribute to the alterations of immune cell subpopulations, including increased type 2 T helper cell and effector memory CD4 T cell, and decreased activated B cell, etc. Given the exploratory nature of this study and the limited validation sample size, these genes should be regarded as candidate biomarkers rather than validated diagnostic or therapeutic targets. Further investigations with larger cohorts and mechanistic studies are warranted to validate their clinical value and biological functions in COPD.