Introduction

Chronic obstructive pulmonary disease (COPD) is a chronic respiratory disorder that primarily manifests as chronic bronchitis and emphysema1 and can progress to a prevalent chronic disease associated with cor pulmonale and respiratory failure.2 The World Health Organization has forecasted that COPD will rank as the third leading cause of death by the year 2030.3 Early identification of chronic respiratory diseases is essential because respiratory conditions remain an important contributor to mortality. Previous study has demonstrated that pulmonary pathologies, including pulmonary edema and pneumonia, constitute a substantial proportion of sudden natural deaths, underscoring the need for improved detection and management of respiratory disorders.4 The majority of patients with advanced COPD die from severe dyspnea, highlighting the importance of early diagnosis and effective treatment. More recently, artificial intelligence has been increasingly explored in the automated diagnostics, clinical decision support, prognostic modelling and precision medicine research of COPD, with great transformative potential.5 While numerous studies have been conducted to diagnose and manage COPD,6 there remain limitations in improving patient survival rates. Consequently, the urgent need arises for the identification of reliable diagnostic markers and promising drug targets to address this pressing issue.

Inflammation is a core feature of COPD, triggering the activation and recruitment of inflammatory cells.7 Natural killer (NK) cells, being innate immune lymphocytes, serve as one of the initial lines of defense against invading pathogens,8 and play a pivotal role in regulating diverse immune responses through the secretion of cytokines and chemokines.9–11 Recent advances have significantly expanded our understanding of NK cell biology in COPD. A recent study has confirmed the essential role of NK cells in COPD pathogenesis and their association with disease progression and exacerbation.12 Mengistu et al have proposed that dysregulated crosstalk between NK cells, dendritic cells, and regulatory T cells significantly contributes to COPD development, highlighting their roles in COPD immunopathology.13 Aberrant NK cell activation has been indicated as a hallmark of advanced COPD, with activated NKp44⁺ NK cells significantly enriched in severe stages.14 While macrophages, monocytes, and neutrophils have been extensively characterized in COPD, NK cells remain relatively understudied despite their critical role as frontline innate immune defenders bridging innate and adaptive immunity. Moreover, TNFAIP3 high NK cells in cancer functionally exhibit immunosuppressive traits.15 It has been indicated that JUNB silence in immune cells could inhibit the IFN-γ expression and secretion from NK cells, leading to reduced STAT1 pathway activation.16 Therefore, elucidating the underlying mechanisms of NK cell function in COPD can offer crucial and novel insights into the therapy for COPD.

Indeed, the significance of NK cells in COPD pathogenesis has been extensively discussed. For example, the NK cell phenotype serves as a crucial indicator for predicting the aggravation of COPD.17 Smoking or exposure to smoke can activate NK cells in the lungs, leading to the accelerated progression of COPD.18 In COPD, dendritic cells enhance the cytotoxic capacity of NK cells via interleukin-15 receptor alpha,19 and the elevated cytotoxicity of NK cells hasten the destruction of lung parenchymal cells,20 additionally, decreased expression of NK cell activating receptors leads to dysfunction of the NK cell population.12 These accumulated evidences suggest that NK cells play a crucial role in the pathogenesis and immune regulation of COPD. Furthermore, several gene expression-based diagnostic models for COPD have been previously developed. For instance, a pyroptosis-related gene signature comprising eight genes has demonstrated diagnostic utility,21 and a macrophage-related gene signature (CLEC5A, FTL, SLC2A3) has been shown to distinguish COPD patients from healthy controls.22 Despite these advances, current diagnostic model studies for COPD face important limitations. Most existing diagnostic signatures have been derived from bulk data analyses without accounting for cellular heterogeneity, and few have focused specifically on NK cell-related genes. Therefore, additional study on NK cells can further enhance our understanding of its pathogenesis as well as develop novel diagnostic model for COPD.

Currently, single-cell RNA sequencing (scRNA-seq) has undergone swift, reliably quantifying transcriptional heterogeneity and enabling comprehensive analyses of both biological systems23 and the immune system.24 In this study, we aimed to address this gap by integrating scRNA-seq to identify NK cell-specific diagnostic signature with bulk RNA-seq validation across multiple independent cohorts (GSE38974, GSE8545, GSE11784 with adequate sample size). After future localized clinical validation, our findings are expected to provide novel insights into the diagnosis and treatment of COPD.

Materials and Methods Subjects

The COPD related bulk datasets GSE38974, GSE8545, GSE11784, and scRNA-seq dataset GSE173896 were downloaded from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The details of all datasets were shown in Table 1 and the available clinical information were summarized in Table S1A–C. The probes in the above datasets were converted to GeneSymbols using the corresponding platform annotation files. Samples containing missing values were excluded prior to subsequent analyses. When multiple probes were mapped to the same GeneSymbol, the average expression value of these probes was adopted.

The ethic approval of our study was exempt from Institutional Review Board of The Third Hospital of Hebei Medical University, based on the item 1 and 2 of Article 32 of the Measures for Ethical Review of Life Science and Medical Research Involving Human Subjects.

Table 1 The Details of Datasets Used in This Study

ScRNA-Seq Data Processing and Analysis

ScRNA-seq data processing was carried out using the R package “Seurat” (version 5.0.1).25,26 Initially, single-cell data were filtered with the following thresholds to remove low-quality cells: total Unique Molecular Identifier (UMI) counts (nCount RNA) > 1000, nFeature RNA (unique molecular identifier, UMI) > 200, nFeature RNA < 5000, and percent.mt <10%. Then, the data were normalized using the “NormalizedData” function. Principal Component Analysis (PCA) was performed on the normalized data using the “RunPCA” function. Unsupervised clustering of the major cell subtypes was accomplished by the “FindClusters” function within “Seurat” package, followed by visualization using tSNE. Next, Human Primary Cell Atlas Database and previously reported NK cell markers (GNLY, GZMB, SAMD3, DOCK2)27 were jointly adopted for NK cell annotation. Pseudobulk processing of NK cells from each sample was performed to alleviate the confounding effects of single-cell sparsity and zero-inflation on differential expression analysis. The AggregateExpression function in the Seurat package was used to construct pseudobulk expression matrices. Briefly, single-cell NK cells were grouped according to sample identity (orig.ident) and phenotypic information (COPD and control groups). The UMI counts of all NK cells originating from the same biological sample were summed to generate expression profiles of pseudobulk samples, and each of which reflected the overall transcriptomic characteristics of the NK cell population in the corresponding individual. The aggregated pseudobulk matrix was finally organized with genes as rows and samples as columns, and the expression values were defined as the total UMI counts of NK cells per sample.

Differential Expression Analysis

To identify the differentially expressed genes (DEGs), differential expression analysis was conducted on the comparison groups using the “limma” R package (version 3.52.4).28 Then the significant DEGs were screened with the cut-off criteria of |log2FC|> 1 and p.adjust <0.05 (p values were adjusted using Benjamini-Hochberg (BH) method).

Functional Enrichment Analysis

The Gene Ontology (GO) terms, Kyoto Encyclopedia of Genes and Genomes(KEGG) pathway enrichment analysis,29 and Gene Set Enrichment Analysis (GSEA) were conducted using the “clusterProfiler” R function package (version 4.7.1.2).30 The p values were determined by Wilcoxon test after Benjamini-Hochberg false discovery rate correction (BH-FDR). A threshold of adjusted p < 0.05 was adopted for identifying significant terms.

Protein-Protein Interaction Analysis

Protein-protein interactions (PPIs) were evaluated using an online STRING database (https://string-db.org/, version 11.0),31 with confidence coefficient >0.4 and subsequently visualized through Cytoscape32 for enhanced comprehension.

Risk Score Model Construction

The LASSO logistic regression analysis can perform variable selection simultaneously with the fitting of generalized linear models.33 After identification of candidate genes, the LASSO logistic regression analysis was performed using the R package “glmnet” (version 4.1.8)34 to screen the hub genes and construct a risk score model. To determine the optimal penalty parameter λ, the five-fold cross-validation was performed to train and validate the LASSO model. Specifically, all samples were randomly divided into five subsets. In each iteration, four subsets were used as the training set and the remaining one as the validation set. After five iterative cycles, the mean cross-validated error (cvm) and corresponding standard error (sd) were calculated. The risk score for each sample was calculated based on the expression of all hub genes:

(1)

where Coefi represented the risk coefficient for each gene, and xi referred to the expression of the gene.

Internal validation was performed to further evaluate the predictive accuracy of the model, and to correct the optimistic bias caused by overfitting. Bootstrap resampling with 200 iterations of sampling with replacement was applied, and the entire procedures of LASSO variable selection and coefficient estimation were repeated in each Bootstrap sample. Furthermore, calibration curves were plotted to assess the consistency between the predicted probabilities and the actual observed outcomes.

Immune Cell Infiltration Analysis

The “CIBERSORT” R package35 was used to calculate the relative proportions of a total of 22 immune cells for each sample in both the high-risk group and low-risk group. A deconvolution algorithm was adopted to characterize the composition of infiltrating immune cells using a predefined signature panel of 547 barcode genes.36 Permutation testing was performed with 1000 permutations to calculate P-values, and quantile normalization (QN) was set as TRUE. The outputs of CIBERSORT encompassed the estimated proportions of each cell subtype, corresponding P-values, correlation coefficients, and root mean square error (RMSE). Only samples with a CIBERSORT deconvolution P-value < 0.05 were retained for subsequent analyses and visualization. The sum of the estimated proportions of all immune cell subsets within each sample equaled 1. In addition, Pearson correlation analysis was performed to assess the association between immune cell infiltration and the risk score. The p < 0.05 was regarded as significant.

Targeted Drug Analysis

The CTD Database (http://ctdbase.org/) was used to assess gene-disease relationships, enabling a deeper understanding of their interconnectedness.

Meanwhile, the DGIdb database (version 4.2.0-sha1 afd9f30b, https://dgidb.genome.wust-l.edu) was utilized to predict targeted drugs for marker genes.

Statistical Analysis

The Wilcoxon test was applied to determine the significance of differences in gene expression and immune cell infiltration between two groups. The “cor” function of R language was performed to analyse the Pearson correlation. Receiver operating characteristic (ROC) curves were generated using the “pROC” R package (version 1.18.4)37 to assess the diagnostic performance of the constructed risk score model and determine its overall diagnostic value.38 All statistical analyses were conducted using R software (version 4.3.1). Statistical significance was considered p <0.05.

Results Single-Cell Clustering and Identification of NK Cell-Related Genes

Using the scRNA-seq data in GSE173896, the cell clustering and annotation was conducted on all cells to identify the cell types in different samples. The results revealed a uniform distribution of cells across all samples, indicating that the data was suitable for subsequent analysis (Figure 1A). Following dimensionality reduction, all cells were categorized into 21 cell clusters (Figure 1B). Concurrently, we examined the distribution of cells across sample types and noted that the cluster composition of COPD samples resembled that of control samples. Nevertheless, we detected disparities in the cell counts in some clusters such as cluster 5, cluster 8, cluster 9, and cluster 11 (Figure 1C).

Figure 1 Identification of NK cell-related genes using single-cell analysis and DEGs analysis in the GSE173896 cohort. (A) A t-SNE plot of cell distribution in seven samples following dimensionality reduction. (B) A t-SNE plot of clustering analysis. (C) The constitution of clusters in COPD and control samples. (D) A t-SNE plot of 11 cell annotation. (E) The proportions of 11 cell types in the COPD and control samples. (F) T-SNE plots of the distribution of 11 cell types in COPD and control samples. (G) Heatmap of DEGs between COPD NK cells and control NK cells. (H) PPI network interaction map of candidate genes. (I) The top 10 KEGG pathway enrichment of NK cell-related marker genes. (J) The GO enrichment of NK cell-related marker genes. The p value was determined by Wilcoxon test, and adjust p value <0.05 was taken as statistically significant.

We identified 11 cell types from the 21 clusters based on marker gene expressions (Figure 1D), including B cells, DC cells, endothelial cells, epithelial cells, fibroblasts, mast cells, monocyte, myofibroblast, neutrophil, NK cells, and T cells. The feature marker expressions of all cell types were displayed in Figure S1. In addition, we analyzed the distribution and proportion of the 11 cell types between COPD patients and controls, and found that the proportion of some cell types was significantly differential, such as endothelial cells, epithelial cells, fibroblasts, and NK cells (Figure 1E and F).

Subsequently, considering the vital role of NK cells in COPD, we conducted a differential expression analysis in the NK cells between COPD and control samples to identify NK cell-related marker genes, a total of 135 DEGs were identified in the COPD NK cells group (Table S2), including 78 up-regulated genes and 57 down-regulated genes compared to the control NK cells group (Figure 1G). Based on these NK cell-related candidate genes, the PPI network included seven genes and nine interactions (minimum required interaction score> 0.15) (Figure 1H).

The KEGG and GO enrichment analyses were performed to obtain the function information of the 135 NK cell-related marker genes. The results showed that 47 KEGG pathways were significantly enriched (p<0.05) (Table S3A), including Th1 and Th2 cell differentiation, and natural killer cell-mediated cytotoxicity, etc (Figure 1I). The GO enrichment results indicated that these marker genes were significantly enriched in multiple important signaling pathways involving inflammatory response and immune regulation, for example T cell differentiation, mononuclear cell differentiation, alpha-beta T cell activation, lymphocyte differentiation, and immune response-activating signaling pathway (Figure 1J and Table S3B).

Identification of Three Candidate NK Cell-Related Genes in COPD

To identify COPD development related genes, the DEGs were further analyzed between COPD and control samples based on bulk data in GSE38974. A total of 624 DEGs were identified in the COPD patients compared to controls (Figure 2A and B), including 311 up-regulated genes and 313 down-regulated genes (Table S4). Next, after taking the intersection of the 624 DEGs and 135 NK cell-related marker genes, three overlapping genes including JUNB, GADD45B, and TNFAIP3 were generated as candidate genes, which were related to COPD development and NK cells in COPD (Figure 2C).

Figure 2 Identification of candidate genes associated with NK cells in COPD. (A) Volcano plots of DEGs between COPD and control samples in the GSE38974 dataset (COPD =23 and control = 9). (B) Heatmap of DEGs between COPD and control samples in the GSE38974 dataset. (C) Venn diagram of the NK cell-related marker genes and DEGs. (D) The top 10 KEGG pathway enrichment of NK cell-related candidate genes. (E and F). The top 10 BP (E) and MF (F) of GO terms of NK cell-related candidate genes. The p value was determined by Wilcoxon test, and adjust p value <0.05 was taken as statistically significant.

We conducted GO and KEGG pathway enrichment analyses to gain insights into the potential functions of these three candidate genes. The KEGG enrichment analysis revealed that they were primarily enriched in 20 pathways like NF-kappa B signaling pathway and TNF signaling pathway (Figure 2D). The Go enrichment results showed that they participated in 229 biological processes (BP) terms such as osteoclast proliferation and regulation of adaptive immune response (Figure 2E), as well as in 12 molecular functions (MF) terms such as K63-linked deubiquitinase activity (Figure 2F). More details of enrichment results were shown in Table S5A and B.

Development of Risk Score Model for COPD Based on NK Cell-Related Genes

To better prevent and indicate the development of COPD, we then constructed a risk model for COPD based on the candidate genes. Using LASSO logistic regression analysis, two genes including JUNB and TNFAIP3 were determined to build a risk score model, based on the minimal criteria of lambda = 2 after five-fold cross-validation (Figure 3A). The risk score model was constructed as: Risk Score = 0.8699024*JUNB+0.7065283*TNFAIP3. The samples in GSE38974 dataset were stratified into high-risk and low-risk groups using the median risk score value as the cutoff. We discovered that the majority of control samples belonged to the low-risk group (Figure 3B). Furthermore, PCA results revealed that high-risk and low-risk samples were relatively separated into two distinct clusters (Figure 3C). These findings strongly suggested that the risk score model was capable of accurately distinguishing COPD samples from controls.

Figure 3 Construction and verification of NK cell-related diagnosis signature in COPD based on the GSE38974 dataset. (A) Selection of the tuning parameter (lambda) in the LASSO model. (B) Distribution of patients based on the risk score. (C) Space distribution of patients based on PCA. (D and E) Box plots of the risk score in COPD and control groups in the GSE38974 dataset (**** p <0.0001) (COPD =23 and control = 9) (D) and GSE11784 dataset (*** p <0.001) (COPD =22 and control = 72) (E). The p values were determined by Wilcoxon rank-sum test. (F and G). ROC curves of the risk score in the GSE38974 dataset (F) and GSE11784 dataset (G). (H). Five-fold cross-validation ROC curves based on GSE38974 and GSE11784 datasets. (I) Heatmap of the expressions of NK cell-related hub genes in high-risk and low-risk groups.

In addition, we have also evaluated the diagnostic performance of the risk model for COPD. Specifically, the GSE38974 dataset was utilized as a training set, and the GSE11784 dataset served as a validation set. Our findings revealed that, in both the training and validation sets, the risk score of COPD patients was significantly higher compared to the controls (Figure 3D and E). This observation suggested that individuals with higher risk scores showed higher risk to develop COPD. Meanwhile, calibration curve and decision curve indicated a favorable reliability of our risk model (Figure S2). ROC analysis showed an AUC value of 0.928 in the training set (Figure 3F) and 0.754 in the validation set (Figure 3G). Furthermore, 5-fold cross-validation ROC curves based on the above two datasets were also analyzed, and both the corresponding average AUC values were more than 0.7, indicating a low risk of overfitting (Figure 3H). These results demonstrated that our risk model exhibited good performance in distinguishing COPD samples from controls. Thus, the subsequent analyses between COPDs vs controls were conducted to characterize the distinct features between high-risk group vs low-risk group samples. Regarding the key genes JUNB and TNFAIP3 in risk model, we noticed that there were significantly differential JUNB, TNFAIP3 expression characteristics between the high-risk group and the low-risk group (Figure 3I).

Significant High Expressions of Two Hub Genes in COPDs Compared to Controls

Furthermore, to validate the expression of two hub genes in COPD, we analyzed their expressions between COPD patients and controls. The results showed that in the GSE38974, GSE8545, and GSE11784 datasets, both of the two hub genes were significantly highly expressed in COPD patients compared to controls (Figure 4A–C).

Figure 4 Validation of the NK cell-related hub genes’ expression in various datasets. (A-C). Violin plot of NK cell-related hub genes’ expression in COPD and controls of GSE38974 (COPD =23 and control = 9) (A), GSE8545 (COPD =18 and control = 18) (B), and GSE11784 (COPD =22 and control = 72) cohorts (C) (p for JUNB <0.05, p for TNFAIP3 p <0.001). All p values were determined by Wilcoxon rank-sum test. * p <0.05, ** p <0.01, *** p <0.001.

Functional Differences Between High and Low Risk Samples

To gain deeper insights into the biological significance of the diagnostic signature, we performed the KEGG and GO enrichment analyses on DEGs between high-risk and low-risk groups using GSE38974. The KEGG enrichment results showed that the DEGs were significantly involved in 35 pathways such as cytokine-cytokine receptor interaction, TNF signaling pathway, HIF-1 signaling pathway, and IL-17 signaling pathway (Figure 5A). The GO terms showed that they were significantly involved in muscle system process, leukocyte migration, etc. (Figure 5B).

Figure 5 Exploration of the biological significance of the risk score model. (A) The top 10 KEGG pathway enrichment of DEGs between high-risk and low-risk groups. (B) The top 10 GO enrichment of DEGs between high-risk and low-risk groups. (C) The activated and suppressed pathways based on GSEA. (D) Several inflammation-related and immune-related pathways were significantly associated with the risk score model. The p values were determined by Wilcoxon test after Benjamini-Hochberg false discovery rate correction (BH-FDR). A threshold of adjusted p < 0.05 was adopted for identifying significant terms.

Furthermore, the GSEA results revealed that various inflammation-related signaling pathways, including the NF−kappa B signaling pathway, TNF signaling pathway, and JAK-STAT signaling pathway, were significantly activated in the high-risk group. Concurrently, certain immune-related pathways, such as B cell receptor signaling pathway and NOD-like receptor signaling pathway, were also notably enriched. On the contrary, the Vascular smooth muscle contraction and ECM-receptor interaction pathways were significantly inhibited (Figure 5C and D). The details of all enrichment analyses were provided in Table S6A–C.

The Divergent Immune Landscapes Between High and Low Risk Groups

Given the crucial role of NK cells in disease immunity, we then investigated the association between our NK cell-related diagnostic signature and immune cell infiltration in COPD. To accomplish this, we employed the CIBERSORT algorithm to analyze the relative abundance of 22 immune cells within the GSE38974 cohort, and observed large differences in immune cell infiltration in each sample (Figure 6A). Furthermore, we found that infiltration proportions of six immune cells, including T cells CD8, NK cells activated, monocytes, macrophages M0, macrophages M2, and mast cells resting, were of significant differences between the high-risk group and low-risk group (Figure 6B). Remarkably, the high-risk group exhibited a significantly increased proportion of monocytes and macrophages M0, whereas the proportion of NK cells activated, T cells CD8, macrophages M2, and mast cells resting were significantly reduced, in comparison to the low-risk group. Furthermore, our data indicated that most of these immune cells also exhibited the similar infiltration tendency in GSE11784 dataset, implying a common infiltration feature of high-risk COPD (Figure 6B).

Figure 6 The landscape of immune cell infiltration in high-risk and low-risk groups. (A) Stacked plot of 22 immune cell types in each sample. (B) The 22 immune cell infiltration between the high-risk and low-risk groups, based on GSE38974 and GSE11784 datasets. (C) Lollipop diagram of the correlation between risk score and immune cell infiltration level; red: statistically significant (p < 0.05). (D) Correlation analysis of 6 vital immune cells and risk score. The p value was determined by Wilcoxon test, and adjust p value <0.05 was taken as statistically significant. * p < 0.05; ** p < 0.01; *** p < 0.001; **** p < 0.0001.

To further validate the association between the differentially infiltrated immune cells and the risk score, we performed correlation analyses. The results demonstrated a significant positive correlation between monocytes (R=0.6, p=0.00028), macrophages M0 (R=0.41, p=0.019), and the risk score (Figure 6C and D). In contrast, NK cells activated (R=−0.38, p=0.03), T cells CD8 (R=−0.51, p=0.0031), macrophages M2 (R=−0.59, p=0.00034), and mast cells resting (R=−0.69, p=0.000013) displayed a noteworthy negative correlation with the risk score (Figure 6C and D).

Prediction of Four Drugs Targeting JUNB and TNFAIP3 in COPD

To further explore the potential clinical value of the NK cell-related hub genes, we have predicted the potential drugs targeting JUNB and TNFAIP3. First, we utilized the CTD database to analyze the inference scores of JUNB and TNFAIP3. The results revealed that among the lung related diseases, both JUNB and TNFAIP3 possessed the strongest correlation with COPD (Figure 7A and B), indicating that our NK cell-related hub genes held promising potential as diagnostic markers for COPD. Subsequently, we delved deeper into exploring potential drugs associated with these two hub genes. Leveraging the DGIdb database, gene-drug interaction results indicated that three drugs were correlated with JUNB and one drug associated with TNFAIP3 (Figure 7C).

Figure 7 Drug prediction of NK cell-related hub genes in COPD. (A and B). Inference score radar map of JUNB (A) and TNFAIP3 (B) in COPD based on the CTD database. (C). NK cell-related hub genes target drug networks predicted by the DGIdb database.

Discussion

COPD is a destructive lung disease ailment marked by chronic inflammation. Over the course of several years, this condition progressively deteriorates, contributing significantly to rising morbidity and mortality rates.1,39 A feature of the chronic inflammatory process in COPD patients is the activation of the innate immune system, with NK cells serving as the frontline defenders against infection, and may be involved in the pathogenesis of COPD.40 Thus, in our present work, we have identified the hub NK cell-related genes and constructed a risk model for COPD, in order to give more insights into their potential associations and to improve early detection of COPD risk.

In this study, we integrated scRNA-seq analysis, bioinformatics analysis, and LASSO logistic regression methods to delve deeply into the transcriptomic data of COPD. Based on COPD related scRNA-seq data, we successfully identified 135 NK cell-related marker genes, which were significantly enriched in immune-related pathways, such as Th1 and Th2 cell differentiation and lymphocyte differentiation. Notably, previous reports indicated that Th1/Th2 cytokines differed in different stages of COPD, and the specific transformation and activation of CD4+ Th1 promoted the immune inflammatory response in COPD patients.41 Additionally, the congenital lymphocyte (ILC) subtype ILC1 was increased in COPD patients.42 Furthermore, both ILC1 and ILC3 were increased in the lungs of patients with severe COPD.43 Hence, these NK cell-related marker genes were probably associated with immune response of COPD.

Going further, after cross-analysis with DEGs basing on bulk data, we screened three COPD NK cell-related candidate genes including JUNB, GADD45B, and TNFAIP3. They participated in inflammatory response pathways, specifically involving the NF-kappa B signaling pathway and TNF signaling pathway. The redox-sensitive transcription factor nuclear factor (NF)-kappa B was recognized as a pivotal player in the regulation of inflammatory responses, playing a pivotal role in both asthma and COPD.44 Additionally, Tumor necrosis factor α (TNF-α), a cytokine renowned for its regulatory role in inflammatory responses, was associated with the origination of several inflammatory conditions.45 Given these connections, we hypothesized that these NK cell-related candidate genes were likely to significantly influence the inflammatory response in COPD.

Utilizing LASSO logistic regression analysis, JUNB and TNFAIP3 were maintained as COPD NK cell-related hub genes and constructed a robust diagnostic signature. While GADD45B is a stress-inducible gene involved in cellular stress responses and inflammatory processes, which may exert biological roles in COPD. Its exclusion suggested that its contribution to diagnostic discrimination was redundant when JUNB and TNFAIP3 were already included. Whether GADD45B plays a role in other aspects of COPD pathogenesis, such as disease progression or exacerbation risk, is deserved to be explored in the future studies. Across multiple datasets, we found that the significantly elevated expression of both hub genes, JUNB and TNFAIP3, in COPD samples compared to controls. Regarding JUNB, while consistently upregulated across all three datasets, it showed variable significance levels, which likely reflected the sample difference and potential impacts of confounders. JUNB, a member of the Jun protein family, played a pivotal role in inflammatory responses, responding to diverse stimuli and regulating the expression of inflammation-related genes. It has been reported that in the inflammatory response of mice during acute infection, JUNB aided T helper 17 cells by regulating AP-1 complex.46 In addition, JUNB was implicated in the differentiation of regulatory T cells.47 Its role in clonal expansion of Th1, Th2, and Th17 cells is also crucial, and JUNB directly inhibited the expression of genes related to apoptosis by enhancing IRF4 DNA binding at the gene locus.48 Notably, JUNB was found to be up-regulated in emphysema lung tissue.49 On the other hand, TNFAIP3 served as a crucial negative feedback regulator, was important to many inflammatory diseases like retinal vasculature,50 and nonalcoholic steatohepatitis.51 During the progression of acute lung injury, it was observed that TNFAIP3 exhibited a significantly over expressed in lung tissue.52 By inhibiting the activation of the NF-κB signaling pathway, TNFAIP3 played a crucial role in inflammatory responses of lung diseases.53 The upregulation of TNFAIP3 observed in our COPD samples presented an interesting biological phenomenon. We speculated that this upregulation likely represented a compensatory anti-inflammatory response aimed at restraining excessive NF-κB activation in the chronic inflammation in COPD. This phenomenon has been observed in other chronic inflammatory conditions, where persistent chronic inflammatory stimuli lead to sustained NF-κB pathway activation that overwhelmed the inhibitory capacity of endogenous regulators such as TNFAIP3.53 In the context of COPD, we thereby hypothesized that the chronic inflammatory milieu drives compensatory upregulation of TNFAIP3, but this response is ultimately insufficient to fully suppress the ongoing inflammatory cascade, resulting in a state of “compensatory exhaustion”. However, the detailed functional role of TNFAIP3 in COPD should be investigated in the future study. While JUNB and TNFAIP3 have been individually implicated in COPD pathogenesis, our study provided a cell-type-specific context for their involvement. These results showed that the up-regulated expression of JUNB and TNFAIP3 in COPDs may activate pulmonary inflammatory processes, and alterations in these genes may be implicated in immunomodulatory mechanisms underlying COPD. Our diagnostic risk model extended the previous knowledge by linking these established inflammatory genes to NK cell biology in COPD.

Based on the above findings, we delved into the relationship between immune cell infiltration and the NK cell-related diagnostic model in COPD. Notably, we discovered a significant positive correlation between the proportion of monocytes and macrophages M0 and the risk score. Monocytes and macrophages M0, renowned as early responders within the immune system, played pivotal roles in initiating the inflammatory response. In monocytes, the inflammatory cascade, primarily driven by NF-κB and TNF signaling, was abnormally activated, subsequently exerting proinflammatory effects on the alveolar epithelial cells of COPD patients.54 Pei et al observed an elevated proportion of monocytes in the peripheral blood of COPD patients.54 Furthermore, studies revealed that in COPD, macrophages exhibit compromised antigen presentation capabilities, reduced cell chemotaxis, and mitochondrial dysfunction, resulting in impaired immune activation among COPD patients.55 In addition, NK cells activated displayed a significant negative correlation with the risk score. NK cells were an innate lymphocyte subpopulation with the ability to influence innate and adaptive immune responses.56 Reduced NK cells may be a result of changes in lymphocyte subsets and immune dysfunction in patients.57,58 The number and function of NK cells were found diminished in severe COVID-19 patients, leading to a decreased elimination of inflammatory cells.59 In summary, these immune cell related evidences in different risk groups suggested that the lungs of COPD patients were enduring a persistent inflammatory response and immune activation process.

Several limitations should be acknowledged in our present study. Although we have integrated multiple single cell and bulk expression datasets from public databases, our analyses remain subject to inherent data limitations, including relatively small sample sizes and heterogeneity across cohorts. Meanwhile, our study is hypothesis-generating and our findings lack experimental validation, thus further wet-lab functional exploration of JUNB and TNFAIP3, along with model validation in larger, more homogeneous cohorts, is warranted in future studies. Moreover, limited by the public data resources, smoking status and other important confounders in COPD have not been included in this work, which should be considered in future validation. Despite the identification of a diagnostic model and potential targets, further large-scale prospective studies are still required before clinical applications, due to lacking generalizable cutoffs and direct mechanistic evidence. Nevertheless, while our study provides a promising NK cell-related risk model, substantial work involving further experimental validation and large-scale clinical studies is warranted before any clinical application.

Conclusions

To summarize, we have established a stable NK cell-related risk score for COPD composed of JUNB and TNFAIP3, which exhibits great potential to distinguish COPD patients from controls. Both hub genes were significantly upregulated in COPD and participated in inflammatory pathways. Our findings provided new insights into the inflammatory and immune mechanisms of COPD and offered potential targets for early diagnosis and therapy for COPD patients. Nevertheless, prospective validation in large, multi-center cohorts and experimental validation in patient-derived samples are essential before further clinical application.

Data Sharing Statement

Data that support the findings of this study are openly available in Gene Expression Omnibus (GEO) at https://www.ncbi.nlm.nih.gov/geo, reference number GSE38974, GSE8545, GSE11784, and GSE173896.

Ethical Approval and Consent to Participate

The ethic approval of our study was exempt from Institutional Review Board of The Third Hospital of Hebei Medical University, based on the item 1 and 2 of Article 32 of the Measures for Ethical Review of Life Science and Medical Research Involving Human Subjects.

Author Contributions

All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.

Funding

This study did not receive support from any organizations.

Disclosure

The authors report no conflict of interest.

References

1. Rabe KF, Watz H. Chronic obstructive pulmonary disease. Lancet. 2017;389(10082):1931–17. doi:10.1016/S0140-6736(17)31222-9

2. Silvestri GA, Young RP. Strange bedfellows: the interaction between COPD and lung cancer in the context of lung cancer screening. Ann Am Thorac Soc. 2020;17(7):810–812. doi:10.1513/AnnalsATS.202005-433ED

3. WHO. The top 10 causes of death. 2024. Available from: https://www.who.int/news-room/fact-sheets/detail/the-top-10-causes-of-death. Accessed August20, 2026.

4. Kanani J. Autopsy analysis of sudden deaths in adults: causes and demographics from a one-year prospective study. Curr Health Sci J. 2025;51(3):343–349. doi:10.12865/CHSJ.51.03.05

5. Chen Z, Tang Z, Ewing RM, Belkhatir Z, Wang Y. Artificial intelligence in respiratory medicine: from diagnosis to treatment and future directions. Chin Med J Pulm Crit Care Med. 2026;4(2):99–116. doi:10.1016/j.pccm.2026.05.005

6. Stolbrink M, Thomson H, Hadfield RM, et al. The availability, cost, and affordability of essential medicines for asthma and COPD in low-income and middle-income countries: a systematic review. Lancet Glob Health. 2022;10(10):e1423–e1442. doi:10.1016/S2214-109X(22)00330-8

7. Caramori G, Casolari P, Barczyk A, Durham AL, Di Stefano A, Adcock I. COPD immunopathology. Semin Immunopathol. 2016;38(4):497–515. doi:10.1007/s00281-016-0561-5

8. Leong JW, Sullivan RP, Fehniger TA. Natural kil