Identification of Oxidative Stress-Associated Biomarkers in Chronic Obstructive Pulmonary Disease: An Integrated Bioinformatics Analysis

Introduction

Chronic obstructive pulmonary disease (COPD) is a heterogeneous ailment characterized by progressively worsening airflow limitation, manifesting clinically through symptoms such as dyspnea, cough, and expectoration.¹ Due to its protracted, recurrent, and evolving nature, COPD has become a prevalent and refractory condition in clinical practice. According to the World Health Organization, the prevalence of COPD is projected to increase significantly over the next four decades, driven by increasing smoking rates in developing nations and aging populations in high-income countries. By 2060, deaths related to COPD and associated disorders are expected to exceed 5.4 million annually.^1,2 This increasing trajectory is intrinsically linked to continuous exposure to COPD risk factors and demographic senescence.³ COPD presents a formidable public health challenge but remains amenable to early preventive and therapeutic interventions. The research elucidates that the etiology of COPD is multifaceted, correlated with smoking, age, utilization of biomass fuels, and genetic predispositions. These risk elements contribute to airway inflammation, oxidative stress (OS) responses, cellular apoptosis, and a protease-antiprotease imbalance, all of which play a critical role in the disease’s pathogenesis. The OS denotes a dynamic disequilibrium between oxidants and antioxidants precipitated by intrinsic or extrinsic stimuli, resulting in damage to DNA, lipids, and proteins, thereby inducing cellular apoptosis and tissue injury.⁴ Chronic inflammation affecting the airways, lung parenchyma, and pulmonary vasculature constitutes characteristic alterations in COPD, where an upsurge in pulmonary reactive oxygen species (ROS) exacerbates pro-inflammatory cytokines, driving the recruitment of inflammatory cells, such as neutrophils and macrophages, augments ROS levels, perpetuating a malignant cycle of OS and chronic pulmonary inflammation. Moreover, OS significantly impacts cellular apoptosis, protease-antiprotease disequilibrium, and immune responses.^5–7 It catalyzes the release of chemokines such as C-C motif chemokine ligand 2, mobilizing dendritic cells and lymphocytes, which partake in immune reactions, thereby promoting autoantibody production and culminating in pulmonary damage. Advancing research has identified several OS-related biomarkers, including glutathione,⁸ protein sulfhydryls,⁹ malondialdehyde,¹⁰ and 8-hydroxydeoxyguanosine.¹¹ While these biomarkers are primarily used in scientific research, they are yet to be widely adopted in clinical practice. Current therapeutic strategies for COPD focus on the administration of glucocorticoids and bronchodilators. The OS-targeted therapies, such as mucolytic agents, including N-acetylcysteine, have demonstrated effectiveness in mitigating acute exacerbations, enhancing patient quality of life, and prolonging survival.¹² Given the critical role of OS in COPD evolution and progression, further exploration into its mechanisms is imperative. This study exploits gene expression data related to COPD from the Gene Expression Omnibus (GEO) database to identify OS-related biomarkers. Using various databases, it performs a comprehensive analysis, establishes and validates a diagnostic model, and predicts relevant microribonucleic acid (miRNA) and targeted therapeutics. These findings aimed to offer potential targets and methods for diagnosing and treating COPD.

Materials and Methods

Data Sources

Gene expression datasets related to COPD were retrieved from the GEO database (https://www.ncbi.nlm.nih.gov/)¹³ using “COPD” as the search term. This search yielded three datasets: GSE130928, GSE11906, and GSE11952. The GSE130928 was designated as the training set, comprising gene expression profiles from alveolar macrophages in bronchoalveolar lavage fluid from 22 patients with COPD and 24 with healthy controls. The GSE11906, serving as the validation set, contained gene expression profiles from the airway epithelium of 33 patients with COPD and 72 with healthy controls. Given that smoking is the primary risk factor for COPD and is closely associated with OS,¹⁴ the GSE11952 was also utilized as a validation set, including gene expression profiles from small airway epithelia of 38 non-smokers (healthy controls) and 45 smokers. The OS-related genes (OSRGs) were acquired from the Human Gene Database (https://www.genecards.org/), totaling 878 genes (Supplementary Table 1). The above data were obtained from public databases, and their download and use were reviewed and approved by the Ethics Committee of the First Affiliated Hospital of Henan University of CM (Opinion No. 2024HL-536).

Construction of Weighted Gene Co-Expression Network and Module Identification

The weighted gene co-expression network analysis (WGCNA) package (version 1.72-5) was used to conduct a WGCNA on the GSE130928 dataset to identify gene modules associated with COPD. Initial steps included clustering of samples to remove outliers and setting the optimal threshold to ensure the WGCNA gene network adhered to a scale-free topology assumption. Based on the similarity of gene expression, genes were categorized into modules, with free genes (genes with low or no change in expression) being discarded. Subsequent association analyses between clinical traits and these modules yielded genesignificance (GS) scores and modulemembership (MM) values for each gene. Finally, modules with a correlation (|Cor|) > 0.4 were selected as key modules for further analysis.

Identification and Functional Enrichment Analysis of Candidate Genes Associated with OS in COPD

Intersections were conducted among differentially expressed genes (DEGs), OS genes, and genes from key modules to identify OSRGs using the Venn diagram package (version 1.7.3). Subsequently, the cluster profiler package (version 4.10.1) was used for gene ontology (GO) biological function enrichment and Kyoto encyclopedia of genes and genomes (KEGG) pathway analyses of the candidate genes. A significance threshold of P < 0.05 was set to filter and visualize enriched biological functions or pathways. Moreover, the OmicCircos package (version 1.32.0) was utilized to depict the genomic locational distribution and expression levels of the candidate genes.

Core Gene Selection and Identification

The STRING database (https://cn.string-db.org/) was applied to establish a protein-protein interaction network, and the network data were imported into the Cytoscape software (version 3.8.2). The cytoHubba plugin was used to analyze the network, determining node degrees to identify core genes. These genes were then validated through receiver operating characteristic (ROC) curve analysis to assess their discriminative power between COPD and control groups within the GSE130928 dataset, with genes achieving an area under the curve (AUC) > 0.7 being designated as core genes. Relative expression levels of these genes in the test and control groups of the GSE130928 dataset were visualized using the ggpubr package (version 0.6.0).

Immune Infiltration and Gene Set Enrichment Analysis (GSEA) of Core Genes

We conducted GSEA enrichment analysis using molecular feature databases To explore the functions of these genes (https://www.gsea-msigdb.org/gsea/msigdb). The downloaded c2. cp. kegg. Hs. symbols. gmt gene set was used as a reference for enrichment analysis of high and low expression samples, with P < 0.05 as the screening condition to determine significantly enriched pathways. The first three pathways with normalized enrichment scores > 0 and < 0 were visualized.

Patients with COPD often have immune dysfunction.¹⁵ We used single-sample GSEA (ssGSEA) to perform immune infiltration analysis on core genes and evaluate the correlation between core genes and 28 types of immune cells.

Construction of Core Gene Regulatory Networks and Drug Prediction

The miRNET database¹⁶ (https://www.mirnet.ca/) was used to predict microRNAs (miRNAs) and transcription factors (TFs) regulating core genes and elucidate the regulatory networks involving core genes and related signal transduction. Furthermore, potential therapeutic drugs for core genes were predicted using the Drug-gene interaction database¹⁷ (https://dgidb.genome.wustl.edu/). Results were visualized using Cytoscape software.

Diagnostic Model Construction and Validation

Constructing a lasso regression diagnostic model using core genes related to OSRG. Gene expression data and gene lists were read, and genes included in the lists were selected from the expression data. The glmnet package (version 4.1-8) was used to apply the least absolute shrinkage and selection operator (LASSO) regression analysis to derive model construction coefficients. To validate the model’s effectiveness, the GSE130928 dataset served as the training set to differentiate between disease and control groups using the model, evaluated using confusion matrix and ROC curve analyses. GSE11906 and GSE11952 datasets were used as validation sets to repeat the analysis and assess the model’s validity.

Results

Identification of DEGs and WGCNA Network Modules

A total of 299 DEGs were identified in the examination of expression data from disease and control groups within the GSE130928 dataset, with 156 upregulated and 143 downregulated genes (Figure 1). To construct the WGCNA network, gene data from all samples were incorporated into a dendrogram, and outliers were excluded. This facilitated the creation of a scale-free network (Supplementary Figure 1), followed by the division of genes into various modules based on expression similarity (Figure 2). A total of 20 modules were identified (Supplementary Figure 2). Subsequent correlation analysis with clinical expressions led to the selection of modules with an absolute correlation coefficient (|Cor|) > 0.4. These modules were labeled with the colors orange, black, dark orange, yellow, sapphire, cyan, and grey (Figure 3).

Figure 1 (a) Differential gene volcano map; (b) Differential gene heatmap.

Figure 2 (a) Gene trees of all samples after removing outliers; (b) Cluster Tree.

Figure 3 Module and Disease Correlation.

Screening and Functional Enrichment Analysis of OS-Related Candidate Genes

By intersecting 299 DEGs with 12,286 module genes and 878 OS genes, 16 OSRGs were identified (Supplementary Table 2). The GO enrichment analysis performed on these OSRGs yielded 734 results, including 645 biological processes (BPs) entries primarily related to OS responses, negative regulation of extrinsic apoptotic signaling pathway, and mediation of apoptotic signaling pathway; a total of 35 cellular component entries mainly involving aggresomes, stress fibers, inclusion bodies, and myosin complexes; and 54 molecular functions (MFs) entries primarily concerning calcium-dependent protein kinase C activity, protein kinase C activity, and calcium-dependent phosphotransferase activity. The KEGG pathway enrichment analysis was used to highlight 38 pathways, involving the AGE-RAGE signaling pathway in diabetic complications, ferroptosis, glutathione metabolism, and the nuclear factor kappa B (NF-κB) signaling pathway (Figure 4). Recent studies have demonstrated that gene loci are closely associated with COPD development and progression.^18,19 Establishing the chromosomal positions of relevant genes helps in understanding whether candidate genes influence COPD through genetic factors. The results indicated that the candidate genes are distributed across various chromosomes (Figure 5), with the most substantial number present on chromosome 6.

Figure 4 (a) Venn diagram of intersected genes; (b) GO enrichment analysis pie chart; (c) KEGG pathway enrichment analysis diagram; (d) GO enrichment analysis bar chart.

Figure 5 (a) Chromosomal positions of candidate genes. Notes: The expression levels of candidate genes in the GSE130928 dataset are represented in the internal circular heatmap. Red represents upregulation, blue represents downregulation, the outer circle of the heatmap represents the control group, and the inner circle represents the disease group; The outermost circle represents chromosomes, and lines from each gene point to their specific chromosomal positions. (b) Manhattan map of candidate genes.

Core Gene Selection and Identification

The 16 OSRGs were identified through the screening. To screen for core genes, a protein interaction network of candidate genes was constructed using the STRING website, and network topology analysis was performed using cytoHubba (Figure 6). Five core genes were identified using the degree value (Supplementary Figure 3): heat shock protein family A (Hsp70) member 1A (HSPA1A; located on chromosome 6), glutamate-cysteine ligase modifier subunit (GCLM; located on chromosome 1), interleukin-1 beta (IL-1β; located on chromosome 2), intercellular adhesion molecule 1 (ICAM1; located on chromosome 19), and glutamate-cysteine ligase catalytic subunit (GCLC; located on chromosome 6).

Figure 6 (a) Network diagram of candidate gene protein interactions; (b) Core gene network diagram.

Subsequent expression analysis of these core genes within the training set indicated elevated expression levels in the COPD group for all five genes (Supplementary Figure 4). The ROC curve analysis yielded AUC values substantiating their diagnostic potential: HSPA1A (AUC = 0.843), GCLM (AUC = 0.841), IL-1β (AUC = 0.860), ICAM1 (AUC = 0.826), and GCLC (AUC = 0.979), all surpassing the threshold of 0.7, thus qualifying for further analytical scrutiny. Co-expression analysis demonstrated a positive correlation among all core genes, with particularly strong correlations between IL-1β and ICAM1, as well as between HSPA1A and GCLM (Figure 7).

Figure 7 (a) ROC curve analysis of core genes; (b) The correlation between core genes and expression levels; (c) IL-1β Correlation with ICAM1 expression; (d) Correlation between HSPA1A and GCLM expression.

GSEA

GSEA was conducted on the core genes to further explore their MFs and associated pathways. All five core genes were enriched in the “cytokine-cytokine receptor interaction” pathway with an upregulated expression profile. Additionally, GCLM, ICAM1, and IL-1β were upregulated in the “chemokine signaling pathway”. The “mitogen-activated protein kinase (MAPK) signaling pathway” indicated upregulation for GCLM, ICAM1, and IL-1β, whereas GCLC was downregulated in this pathway. In the “Parkinson’s disease pathway”, GCLC, GCLM, and IL-1β were all downregulated. Furthermore, low expression levels of HSPA1A, ICAM1, and IL-1β were associated with the “ribosome” pathway (Figure 8).

Figure 8 Gene set enrichment analysis (a) HSPA1A; (b) GCLM; (c) IL-1β; (d) ICAM1; (e) GCLC.

Immune Cell Infiltration Scoring and Core Gene Correlation Analysis

A delineation of associations between pivotal genes and diverse immune cellular phenotypes was established using immunocyte infiltration scores. The HSPA1A was used to manifest significant positive associations with adipocyte-like mast cells, nascent B cells, plasmacytoid dendritic cells, central memory CD4+ T cells, activated CD4+ T cells, type 2 helper T cells, and gamma-delta T cells while exhibiting a notable negative correlation with activated B cells. The GCLM displayed pronounced positive correlations with adipocyte-like mast cells, activated CD4+ T cells, and plasmacytoid dendritic cells and was inversely related to regulatory T cells (Tregs), T follicular helper cells, activated B cells, macrophages, and myeloid-derived suppressor cells. The IL-1β was positively correlated with adipocyte-like mast cells, cytotoxic natural killer T cells, activated CD4+ T cells, natural killer cells, nascent B cells, plasmacytoid dendritic cells, and type 17 helper T cells (Th17), and negatively with regulatory T cells. The ICAM1 correlated positively with adipocyte-like mast cells, monocytes, plasmacytoid dendritic cells, activated CD4+ T cells, natural killer cells, activated dendritic cells, central memory CD4+ T cells, and Th17 cells. It demonstrated negative correlations with regulatory T cells, activated B cells, and immature dendritic cells (Figure 9). The GCLC did not exhibit notable correlations with any immune cells. These findings underscored the significant relationships between core genes and a spectrum of T cells, B cells, dendritic cells, and mast cells, corroborating previous research on their roles in COPD pathogenesis.²⁰

Figure 9 Correlation between immune cell infiltration score and core genes (a) HSPA1A; (b) GCLM; (c) IL-1β; (d) ICAM1. Note: The highlighted red part on the vertical axis in the figure represents the correlation between genes and this type of immune cell (P<0.05).

Core Gene Predictive Drug, miRNAs, and TFs Analysis

The Dgidb database was used to facilitate the prognostication of 33 pharmacological agents targeting these core genes (Supplementary Table 3). Specifically, 20 agents were predicted for IL-1β, including erythromycin, aspirin, hydrocortisone, ofloxacin, cephalexin, and gevokizumab; 5 agents including pegolimab and BI-505, were predicted for ICAM1; 4 agents including cisplatin and sulfamethoxazole were foreseen for GCLC; and predictions for GCLM and HSPA1A involved 3 and 1 agents respectively. The miRNet database was employed to predict 33 miRNAs (Supplementary Table 4) and 37 TFs (Supplementary Table 5) impacting these core genes, excluding GCLM (Figure 10). For instance, miR-335-5p was identified as a regulator for HSPA1A, ICAM1, and GCLC; TFs such as NF-κB1 and RelA were shared regulators for GCLC, ICAM1, and IL-1β.

Figure 10 (a) Drug prediction affecting core genes; (b) Prediction of miRNAs regulating core genes; (c) Predicting transcription factors that regulate core genes.

Diagnostic Model Establishment and Validation for Core Genes

A diagnostic model predicated on these five core genes was constructed using LASSO regression and ROC analysis, with the trajectory of regression coefficients depicted. Subsequent validation through the ROC curve evidenced an AUC of 0.981, affirming the model’s diagnostic precision. Besides, a confusion matrix demonstrated an accuracy of 0.93, precision of 0.96, recall rate of 0.92, and specificity of 0.95 (Figure 11). To corroborate this model, GSE11906 served as an external validation set for disease, indicating an AUC of 0.949, with the confusion matrix reporting accuracy of 0.94, a precision of 0.92, a recall rate of 1.00, and a specificity of 0.81. The GSE11952, used as an external validation set for risk factors, indicated an AUC of 0.906, with the confusion matrix revealing an accuracy of 0.87, precision of 0.85, recall rate of 0.87, and a specificity of 0.87 (Figure 12). These metrics underscore the model’s general applicability and potential as a pre-diagnostic tool for diseases.

Figure 11 Establishment and validation of LASSO regression model; (a) 10 fold cross validation chart; (b) Regression coefficient path diagram; (c) Training set ROC curve; (d) Training set confusion matrix.

Figure 12 (a) ROC curve of disease validation set; (b) Disease validation set confusion matrix; (c) Risk factor validation set ROC curve; (d) Risk factor validation set confusion matrix.

Discussion

COPD is a common respiratory system disease, with OS as its primary pathogenic mechanism. Currently, clinical interventions targeting OS are limited, and there is a scarcity of relevant diagnostic biomarkers. Early diagnosis or screening of COPD remains a challenging issue in clinical settings. With the recent escalation in environmental pollution, including smoke from cigarettes and cooking oils, and PM2.5, there has been an increase in both exogenous ROS caused by these pollutants and endogenous ROS produced by lung inflammation and structural cells. This increase results in damage to the airway epithelial cells, induces the proliferation of smooth muscle cells, leads to a decline in lung function, and makes the OS response increasingly significant in COPD development and progression.²¹ As a result, selecting biomarkers related to OS and the development of drugs targeting these markers might become hotspots in future COPD diagnosis and treatment research. In this study, we identified 16 biomarkers related to OS through the screening of DEGs in the GSE130928 dataset. The GO enrichment analysis indicated that these OSRGs are mainly involved in BPs such as OS responses, regulation of apoptotic signaling pathways, and protein kinase activity. Previous research has revealed that an enhanced OS response intensifies the process of cellular apoptosis.²² The KEGG pathway enrichment analysis was used to demonstrate that OSRGs are primarily involved in pathways, including the AGE-RAGE signaling pathway related to complications in diabetes, ferroptosis, glutathione metabolism, and the NF-κB signaling pathway. These pathways are highly associated with OS, such as the activation of the NF-κB pathway, which induces the expression of various pro-inflammatory cytokines such as tumor necrosis factor-alpha, IL-6, and IL-1β, thereby promoting apoptosis and the OS process.²³ To further analyze the functions of OSRGs, we determined their chromosomal positions and observed that OSRGs are related to multiple chromosomes. Relevant studies indicate that the X chromosome is associated with COPD-related phenotypes, and a variant near TMSB4X, rs5979771, reveals genome-wide significance with lung function.¹⁹ Our results indicate that most genes are located on chromosome 6, and previous studies have also confirmed that certain single nucleotide polymorphisms (SNPs) on chromosome 6 are consistently associated with early susceptibility to COPD.²⁴ These findings suggest that these chromosomes may be significantly related to COPD, and identifying COPD-related SNPs also has certain clinical significance.

Through further selection, we identified five core OSRGs: HSPA1A, GCLM, IL-1β, ICAM1, and GCLC. Based on these genes, we explored their potential MFs through GESA, immune infiltration scoring, miRNAs, and TF predictions while developing a diagnostic model predicated on these core OSRGs, characterized by high sensitivity and selectivity. HSPA1A, GCLM, and GCLC consistently exhibited elevated expression levels in COPD samples. The HSPs form a superfamily whose proteins may increase in response to cellular stress associated with pollutants. Previously reported SNPs within HSP genes are linked to the risk and severity of COPD, and intracellular HSP levels may vary with different external exposures.²⁵ Moreover, a study involving coal miners indicated that elevated plasma levels of HSPA1A could be associated with an increased risk of COPD among these workers.²⁶ GCLM and GCLC are part of the glutathione reductase system, intimately involved with the synthesis of glutathione in the body. Research has indicated that polymorphisms at the GCLM gene locus are related to COPD susceptibility, and GCLC expression is significantly upregulated in both acute exacerbation of COPD and stable patients with COPD, a phenomenon possibly linked to increased methylation of GCLC.^27,28 IL-1β is a pro-inflammatory cytokine, and ICAM1 is a critical inflammatory mediator. Over-secretion of IL-1β can induce ICAM1, thereby exacerbating inflammation.²⁹ Studies indicate that during acute exacerbations of COPD, serum levels of IL-1β and IL-17 are significantly higher than in stable COPD or control groups, correlating positively with serum C-reactive protein levels, neutrophil percentages, and smoking status.³⁰ These findings suggest that increased levels of IL-1β and ICAM1 are linked to pulmonary oxidative damage and are positively correlated with inflammation levels. Another study confirmed that increased expression of IL-1β is closely associated with smoking in both smokers and non-smokers from healthy populations.³¹

Further exploration of the core genes revealed that they are enriched in the chemokine and the MAPK signaling pathway, which is persistent in previous studies.^32,33 The MAPK is critical in the physiological and pathological development of COPD by activating key TFs and inducing the expression of cytokines and chemokines. Immune cell infiltration scoring has elucidated the correlation of core genes with immune cells, including mast cells, CD4+ T cells, and dendritic cells, which aligns closely with previous research.^12,34–36 This provides a reference for subsequent treatment strategies and mechanistic explorations for COPD. Predictions of drugs, miRNAs, and TFs indicated various antibiotics, antibodies, and chemotherapy agents that could intervene with the core genes. However, no related drugs were predicted for GCLM, representing a potential direction for future research. Currently, using antioxidants in clinical settings primarily revolves around various expectorants, which, although they can improve and delay acute exacerbations of COPD, whether this effect is due to the mitigation of OS still requires further investigation.^37,38 Consequently, it is necessary to develop new antioxidants or targeted drugs to provide new options for the clinical treatment of COPD. We predicted 33 miRNAs related to COPD, such as miR-203a-3p, which is highly expressed in patients with COPD and smokers and is associated with basal cell proliferation.³⁹ The miR-221-3p can alleviate cell apoptosis and inflammatory responses in COPD.⁴⁰ Among TFs and NF-κB, a dimeric TF involved in inflammation, immune response, and cell proliferation belongs to the Rel protein family. Another member of this family, RelA, was also predicted, and both are implicated in regulating genes involved in the COPD process, a mechanism that has been confirmed in past studies.^41,42 The diagnostic model constructed based on the core genes demonstrated high accuracy and sensitivity in distinguishing patients from controls in both training and validation datasets. Additionally, attempts to use this model to differentiate between smokers and non-smokers among healthy individuals demonstrated high accuracy, which might be linked to the OS induced by smoking. Therefore, this model may have clinical value for diagnosing COPD and the potential for early screening of the disease.

Conclusion

In this study, we utilized bioinformatics tools to identify five genes highly associated with OS during COPD progression. We comprehensively unveiled the relationship between these genes and the disease, along with their potential for clinical application translation through subsequent in-depth analysis of MFs, BPs, immune infiltration, drug predictions, miRNAs, and TFs, as well as the construction of a diagnostic model. However, further experimental validation is still required to clarify the OS functions related to COPD, including the mechanisms of many upstream regulatory miRNAs that remain undefined. Moreover, the function of this diagnostic model to differentiate between smokers and non-smokers still needs further confirmation and validation. Overall, while the steps in this study were relatively comprehensive, there are shortcomings, such as the small sample size of the dataset and the limited number of OSRGs identified. Nonetheless, the high sensitivity and accuracy of the model provide new references for COPD clinical diagnosis and offer direction for subsequent research and the development of related pharmaceuticals.