Single Cell Transcriptomics Genomics Based on Machine Learning Algorithm: Constructing and Validating Neutrophil Extracellular Trap Gene Model in COPD

Introduction

Chronic obstructive pulmonary disease (COPD) is now recognized as a complex, multicomponent disease characterized by chronic systemic inflammation and is currently the third leading cause of death worldwide.¹ The mortality rate for male COPD patients is higher, with an increasing trend in mortality rates among those over 45 years of age.² In developing countries, mortality rates are also increasing in tandem with the rise in smoking rates. In China, tobacco-related deaths account for 12% of all deaths, and predictions suggest that this proportion could reach 33% by 2030.³ Although the hazardous role of smoking in the development of COPD is well established, nearly half of COPD patients are non-smokers, particularly women exposed to biomass smoke in poorly ventilated homes.⁴

The high mortality rate associated with COPD is linked to its progressive, irreversible airway obstruction and complex comorbidities.^5,6 Each bacterial or respiratory viral infection experienced by this population exponentially increases the risk of adverse outcomes.⁶ Currently, identified endogenous factors associated with COPD include alpha-1 antitrypsin (AAT) deficiency and geNETic mutations in the telomerase gene TERT, which can lead to early-onset emphysema and are also seen in idiopathic pulmonary fibrosis.^7–9

The role of inflammatory mechanisms in COPD has been progressively investigated, characterized by increased and activated macrophages and neutrophils.¹⁰ Sputum neutrophilia is a significant feature of COPD, and this neutrophilic inflammation is induced by cigarette smoke, bacteria, viruses, and oxidative stress.¹⁰ The increase and local chemotactic activation of neutrophils are associated with resistance to inhaled corticosteroids, which becomes a critical factor in exacerbating COPD mortality and progression.¹¹ Although CXCR2 antagonists (such as navarixin) can block the chemotactic effects of CXCL8 and related chemokines to reduce the number of neutrophils in the sputum of COPD patients, they do not improve already deteriorated lung function.¹²

In recent years, the formation of neutrophil extracellular traps (NETs) has gained increasing attention. Neutrophils, when recruited, have a certain probability of undergoing programmed cell death, leading to the rupture of cell membranes and the release of web-like structures containing extracellular double-stranded DNA (dsDNA) combined with histones, myeloperoxidase, and neutrophil elastase.¹³ In respiratory viral infections, NETs formation is considered one of the beneficial protective mechanisms.¹⁴ The formation of NETs is advantageous for the entrapment and containment of pathogens, and it has been reported that NETs exhibit antiviral properties against poxviruses.¹⁵ Conversely, persistent NETs formation is considered a driver of immune pathology, most extensively demonstrated in COVID-19, where systemic and airway NET accumulation induces high inflammation in both acute and chronic disease stages.^16,17 NETs, loaded with highly basic histones and degradative enzymes, possess cytotoxic effects. When they are formed in excess or not cleared by the still poorly understood mechanisms, they can directly lead to host cell death and chronic tissue damage.¹⁸ The tissue damage caused by local overload of NETs is closely related to pulmonary fibrosis, acute respiratory distress syndrome (ARDS), and the progression of asthma.^18,19 Studies have found that the sputum of exacerbated COPD patients contains a large amount of NET.²⁰ Importantly, NETs formation is not limited to exacerbated COPD; it also exists in stable COPD patients and correlates with the severity of airflow limitation.²⁰ Pharmacologists have been committed to manipulating NETs formation and reducing the pathogenic effects of NETs by administering DNase to degrade the NETs scaffold.²¹ In recent years, with the development of genomics, a series of gene signatures related to NETs formation have been gradually identified in various diseases and have been well applied in predicting ischemia-reperfusion injury after lung transplantation, sepsis, and thrombotic diseases.^22–25

Hematological and immunological data (HID) are crucial for the diagnosis, prognosis, and monitoring of numerous diseases.²⁶ As a result, HID data have been widely employed in artificial intelligence research to diagnose and predict the course of diseases.^27–29 In recent years, various clinical studies have demonstrated that HID is an effective and cost-efficient alternative for the early diagnosis and prognosis of coronavirus disease 2019 (COVID-19) and other diseases.^29–31 For instance, Huyut et al successfully developed artificial intelligence models using HID data for the diagnosis and prognosis of COVID-19.^32–34 Yan et al effectively applied the Gradient Boosting Decision Tree method to HID data for the early diagnosis of multiple myeloma.³⁵ Soerensen et al developed an artificial intelligence model based on HID to predict cancer risk in primary care.³⁶ Wu et al achieved 95.7% accuracy in diagnosing lung cancer using the Random Forest algorithm based on 19 routine blood test parameters.³⁷ Haider and others used artificial neural NET works to analyze HID data for early differentiation among leukemias.³⁸

Compared with existing studies, this study innovatively identified NET-related genes associated with poor COPD prognosis using three machine learning methods combined with single cell sequencing dataset and multiple large sample lung tissue transcriptome sequencing dataset. We also constructed prognostic prediction models suitable for both smoking and non-smoking patients and preliminarily confirmed their good efficacy. This research will provide a basis for understanding the pathogenesis of COPD driven by gene-mediated NETs formation and for targeting NETs to improve the prognosis of COPD patients.

Materials and Methods

Data Acquisition and Preprocessing

We downloaded single-cell RNA sequencing (scRNA-seq) data of COPD lung tissue (GSE173896) and transcriptome sequencing data of COPD (GSE37768, GSE54837, GSE57148), and ARDS (GSE76293) from the Gene Expression Omnibus (GEO) database (www.ncbi.nlm.nih.gov/geo/). Ultimately, we obtained scRNA-seq profiles from five COPD patient lung tissues and three age-matched non-COPD lung tissues as controls from GSE173896. The remaining transcriptome data were used to subsequently identify differentially expressed genes. The whole research process is shown in Supplemental Figure 1. In the entire workflow, first, the neutrophil signature genes in COPD tissues were successfully identified in comparison with normal tissues. Meanwhile, the following intergroup comparisons were conducted on the COPD-RNA-Seq data: high GOLD stage vs low GOLD stage, smokers vs non-smokers, and ARDS occurrence vs non-occurrence groups, to identify DEGs for subsequent machine learning. Subsequently, the preliminary validation efficacy was assessed using alveolar lavage fluid samples from COPD patients.

Acquisition of COPD Neutrophil Differential Gene Set

In order to obtain the genes associated with the neutrophil extranuclear trap, we first need to perform cluster analysis of single-cell sequencing samples of COPD and successfully annotate the neutrophil subsets. First, we performed quality control on the downloaded GSE173896 dataset to ensure the accuracy and reliability of the data. This included checking sequencing depth, excluding low-quality cells, and removing cells with a high proportion of mitochondrial genes (Supplementary Figure 2). We created Seurat objects based on the following filtering criteria: genes expressed in at least 3 cells and cells with at least 200 genes expressed. At the same time, samples with a total UMI number less than 1000 and samples with a gene number less than 500 were filtered out. The first 2000 highly variable genes were then selected for subsequent analysis. We utilized the RunPCA function from the Seurat package to perform principal component analysis (PCA) on the preprocessed data, selected high variable genes (HVGs) for dimensionality reduction using PC4, and employed the Run-TSNE function for t-SNE clustering analysis. The results of t-SNE cell subsets clustering in all samples of COPD group and normal control group are shown in Figure 1A and B. Clustered cell populations were annotated using Cell Marker 2.0 to identify components.³⁹ Differential analysis with control group data identified 165 highly expressed neutrophil differential genes.

Figure 1 Single cell sequencing analysis of COPD and Control samples. (A) Preliminary results of tSNE clustering using high-mutation genes from PC4 in dimensionality-reduced data obtained from five COPD samples and three control samples. (B) Twenty-four cellular clusters identified through clustering at the same resolution. (C) Volcano plot of differential genes within the neutrophil cluster, with blue representing downregulated differential genes and red representing upregulated differential genes. (D) tSNE plot after cell cluster annotation. (E) Stacked plot of annotated cellular components, with different colors representing distinct cell clusters.

Identification of Differentially Expressed Genes (DEGs) Between Different Subgroups of COPD

Having successfully identified a range of neutrophil subpopulation specific genes from single-cell cluster analysis, we then need to identify which genes are specific to COPD disease, which genes are associated with poor disease prognosis, and even which genes may be associated with smoking as a risk factor. After obtaining gene mRNA expression matrices and clinical information for 18 COPD samples, 11 smoking healthy controls, and 9 non-smoking healthy controls from GSE37768; 68 GOLD stage 3/4 (stage high), 68 GOLD <3 (stage low), 84 smoker controls, and 6 non-smoker controls from GSE54837; and 12 ARDS samples and 12 normal controls from GSE76293, we preprocessed the data, including gene matching, removal of blank information, and standardization. We used the R package “limma” for intergroup differential analysis and “ggplot2” for visualizing DEGs, setting the criteria to |log2 fold change (FC)| > 1.0 and adjusted p-value < 0.05. Additionally, we employed the R package “cluster Profiler” for KEGG and GO enrichment analysis of differential genes. Furthermore, we used the MCP counter estimate function to analyze the infiltration of eight common immune cells in the GSE37768, GSE54837, and GSE76293 datasets.

RF-LASSO-SVM RFE Machine Learning for Feature Gene Selection and Model Construction

Machine learning can better identify the weight of different factors in different disease states. Therefore, we carry out machine learning on the differential genes obtained in the previous step to identify which genes are the most critical. We intersected the COPD smoking DEGs and COPD nonsmoking DEGs obtained above with COPD stage high DEGs, ARDS DEGs, and scRNA-NDEGs. Finally, we extracted 23 smoking COPD-NDEGs and 33 nonsmoking COPD-NDEGs and included them in the random forest algorithm, LASSO regression analysis, and SVM RFE machine learning. The three machine learning methods were completed using the R packages “randomForest”, “glmnet”, and “e1071”, respectively. All variables were used for SVM model construction, with AvgRank sorted by the average ranking of 10-fold cross-validation. Ultimately, we obtained 3 COPD smoking NDEGs (CD63, RNASE2, and ERAP2) and 4 COPD nonsmoking NDEGs (GRIPAP1, NHS, EGFLAM, and GLUL) for the construction of models A and B through intersection.

Single-Gene Efficacy Assessment and Model Efficacy Validation

We used the GSE37768 dataset to perform group expression analysis and obtain ROC diagnostic dependency curves for the 7 feature genes selected by machine learning to assess the efficacy of single genes in identifying COPD. Subsequently, we used the R package “clusterProfiler” for GSEA analysis of highly expressed differential genes. Additionally, we studied the protein molecular interaction patterns of the two groups of genes using protein-protein interaction (PPI) NETworks and STRING database-based protein interaction NETwork analysis. Furthermore, we constructed gene diagnostic nomograms for the two groups of NDEGs. To validate the efficacy of the machine learning-built models A-B, we used GSE37768 as the internal training set and GSE54837 and GSE57148 as external validation sets for risk score grouping validation and construction of corresponding ROC curves.

Distribution and Protein Structure Prediction of RNASE2 and NHS

We investigated the distribution of RNASE2 and NHS genes through the website Human Protein Atlas (proteinatlas.org), including mRNA distribution in different human organs and various immune cell subsets, as well as subcellular localization prediction. Finally, we constructed protein structure models for RNASE2 and NHS using Alphafold 2.0.

qRT-PCR and ELISA Experiments on COPD Bronchoalveolar Lavage Fluid Samples

We collected bronchoalveolar lavage fluid from ten COPD patients for the detection of RNASE2 and NHS genes and MPO protein ELISA. Among them, five patients requiring ventilator support were classified as stage high, and the other five patients not requiring ventilator support were classified as stage low. The sequence of primers associated with each gene for qPCR is provided in Supplemental Table 1. mRNA was isolated using a commercial Total RNA Purification Kit from NORGEN, and its concentration was determined by a Thermo Fisher Scientific spectrophotometer. The mRNA was then converted to cDNA using the High-Capacity RNA-to-cDNA Kit from Applied Biosystems. Quantitative PCR (qPCR) was performed with SYBR green PCR master mix, also from Applied Biosystems, on a Bio-Rad CFX Manager system. Finally, Δct values were used as the relative expression of RNASE2 and NHS genes compared to the GAPDH reference gene. A commercial MPO ELISA kit was used, which includes a pre-coated microplate with specific anti-MPO antibodies, enzyme-labeled secondary antibodies, washing buffer, substrate solution, and stopping solution. According to the kit instructions, the supernatant samples obtained after centrifugation of bronchoalveolar lavage fluid were appropriately diluted, and standard curves were constructed using standard samples at different gradient concentrations. Then, experimental samples were added, incubated, washed, and developed according to the instructions, and finally, the absorbance (OD values) of each well was measured at 450 nm using an enzyme-linked immunosorbent assay reader. The obtained sample OD values were converted into relative concentrations using the standard curve.

Statistical Methods

Data analysis and visualization were conducted using R software (version 4.3.3) along with GraphPad Prism 8. For the group difference analysis of RT-PCR and ELISA data, we employed both the two-tailed unpaired Student’s t-test and the Wilcoxon rank-sum test. For the differential analysis of GEO transcriptome data, we utilized the limma-voom method to assess differences between groups. A statistically significant result was considered when the P-value was less than 0.05.

Results

Successful Identification of COPD Neutrophil Differential Genes via scRNA-Seq

We analyzed the scRNA-seq profiles of lung tissues from five COPD patients in the GSE173896 dataset, with three age-matched non-COPD lung tissues as controls. As shown in Supplementary Figure 2, after rigorous quality control analysis, data standardization and scaling, and selection and filtering of highly variable feature-expressing cells, four principal component analysis modules were successfully constructed (Supplementary Figure 2F), with the top 20 genes of each PCA displayed in Supplementary Figure 2G. In Figure 1A and B, following tSNE nonlinear dimensionality reduction of PC4, 24 cellular clusters were preliminarily identified. Neutrophil populations in both COPD and control samples were successfully identified after cellular annotation, yielding 165 overexpressed neutrophil differential genes (scRNA-NDEGs) (Figure 1C). The composition and proportion distribution of cells are shown in Figure 1D and E, with a significantly higher neutrophil ratio in COPD samples compared to controls.

Acquisition of COPD and ARDS Differential Genes from Large-Scale Transcriptome Datasets

To obtain differential genes between COPD and smoking and non-smoking controls, we performed group differential analysis on the mRNA data of 18 COPD samples, 11 smoking healthy controls, and 9 non-smoking healthy controls from the GSE37768 dataset, yielding 178 COPD-Smoking-DEGs and 511 COPD-Non-smoking-DEGs. To assess differences among COPD patients of different GOLD stages, we conducted group differential analysis on 68 GOLD stage 3/4 (stage high), 68 GOLD <3 (stage low), 84 smoker, and 6 non-smoker controls from the GSE54837 dataset. Finally, 708 COPD-stage high DEGs, 2076 COPD stage high-smoking-DGEs, and 904 COPD stage high-Nonsmoking-DEGs, were identified. To further analyze gene differences in lung tissue samples as COPD progresses to ARDS, we performed differential analysis on 12 ARDS samples and 12 normal control samples from the GSE76293 dataset. We successfully identified 3320 DEGs in the transcriptome sets of ARDS and normal control samples. We selected the top 30 differential genes from the above subgroup analysis and displayed them in the form of heat map (Supplementary Figure 3A–F). Meanwhile, the results of gene expression difference analysis for each subgroup were displayed in the form of volcano map (Supplementary Figure 4A–F).

Candidate Gene Selection for Smokers and Non-Smokers and KEGG and GO Enrichment

To explore the impact of smoking and non-smoking on neutrophil proliferation, activation, and disease progression in COPD patients, we intersected the multiple transcriptome differential genes obtained with scRNA-NDEGs. Results in Supplementary Figure 5G and H show that we identified 23 smoking COPD-NDEGs and 33 non-smoking COPD-NDEGs for subsequent machine learning. Additionally, we performed KEGG and GO enrichment analyses on the differential genes obtained from the transcriptome samples, finding that genes highly expressed in COPD stage high were enriched in membrane potential regulation, passive transmembrane transporter activity, and ion channel activation pathways (Supplementary Figure 5E and F). COPD stage high differential genes are closely related to the formation of cell membranes and ion channel complexes (Supplementary Figure 5E and F). Compared to the non-smoking group, COPD-Smoking-DEGs were mainly enriched in DNA and cell signal transduction pathways, phosphatidylcholine metabolism pathways, and DNA damage repair checkpoint pathways (Supplementary Figure 5A and B). In contrast, COPD-NonSmoking-DEGs were primarily enriched in pathways related to endothelial cell proliferation and migration (Supplementary Figure 5C and D).

MCPcounter Immune Infiltration Analysis

To assess differences in immune cell infiltration among different COPD subgroup samples, we used the “MCPcounter” R package to quantify and score the distribution differences of eight immune cells in the GSE37768, GSE54837, and GSE76293 datasets. Results in Supplementary Figure 6A and B show that the COPD group had higher infiltration of T cells, monocytes, and neutrophils, with the most significant difference in neutrophil infiltration. Compared to the non-smoking group, smoking group samples showed a trend of increased neutrophil and monocyte infiltration. The COPD stage high group had significantly higher neutrophil infiltration than the COPD stage low and other control groups. Interestingly, higher neutrophil infiltration scores were also observed in ARDS samples (Supplementary Figure 6C).

Machine Learning Model Construction Using RandomForest, LASSO, and SVM-REF

To evaluate the disease-related effects of the 23 smoking COPD-NDEGs and 33 nonsmoking COPD-NDEGs obtained in previous steps and to construct an efficient model, we employed three machine learning methods: RandomForest, LASSO, and SVM-REF, to jointly filter key feature genes. For the smoking group, distinct outcomes were observed: The SVM-REF model prioritized 5 genes to optimize performance metrics (AvgRank < 10), as shown in Figure 2A and B. LASSO regression analysis, validated through 10-fold cross-validation, proposed 3 genes (CD63, RNASE2, and ERAP2) for Model A construction (Figure 2C and D). Gene importance rankings derived from the RandomForest algorithm are presented in Figure 2E and F. The intersection of gene selections across methods was visualized in Figure 2G, ultimately defining the final candidate genes for downstream modeling. For the nonsmoking group, the SVM-REF model identified 7 genes required to achieve maximum accuracy and minimum error rate under the condition of AvgRank < 10, as illustrated in Figure 3A and B. Following 10-fold cross-validation, LASSO regression further refined the selection, recommending 4 genes (GRIPAP1, NHS, EGFLAM, and GLUL) for constructing Model B (Figure 3C and D). The importance rankings of these genes, evaluated by the RandomForest algorithm, are displayed in Figure 3E and F. A Venn diagram analysis integrating results from the three methods (Figure 3G) highlighted candidate genes through intersection screening, which were subsequently used to build gene feature models. ALL the candidate genes after intersection screening were constructed by following gene feature models:

Figure 2 Machine Learning Analysis and Model Construction with Smoking COPD-NDEGs. (A and B) display the results obtained using the SVM-RFE algorithm, illustrating the number of feature genes recommended for achieving the highest accuracy and the lowest error rate, respectively. (C and D) present the outcomes of incorporating feature genes into LASSO regression analysis. The y-axis in (C) represents the coefficients assigned to each feature gene, while the dashed line in (D) indicates the optimal number of factors recommended. (E and F) showcase the results derived from stepwise random forest analysis. The lollipop chart in E is a sorted display based on the importance coefficients assigned. (G) is a Venn diagram representing the intersection of genes recommended by three machine learning methods. The machine learning above uses the smoking COPD-NDEGs obtained in the previous step.

Figure 3 Machine Learning Analysis and Model Construction with nonsmoking COPD-NDEGs. (A and B) display the results obtained using the SVM-RFE algorithm, illustrating the number of feature genes recommended for achieving the highest accuracy and the lowest error rate, respectively. (C and D) present the outcomes of incorporating feature genes into LASSO regression analysis. The y-axis in (C) represents the coefficients assigned to each feature gene, while the dashed line in (D) indicates the optimal number of factors recommended. (E and F) showcase the results derived from stepwise random forest analysis. The lollipop chart in E is a sorted display based on the importance coefficients assigned. (G) is a Venn diagram representing the intersection of genes recommended by three machine learning methods. The machine learning above uses the nonsmoking COPD-NDEGs obtained in the previous step.

Model A (smoker): risk score = (0.126 * expression of CD63) + (0.535 * expression of RNASE2) + (0.019 * expression of ERAP2)

Model B (nonsmoker): risk score = (0.353 * expression of GRIPAP1) + (1.497 * expression of NHS) + (0.0726 * expression of EGFLAM) + (0.627 * expression of GLUL)

Single-Gene Evaluation

To assess the feature genes selected by machine learning, we used the GSE37768 dataset as a training set for single-gene efficacy validation. Evaluation results of the 3 candidate genes for Model A (Figure 4A–C) showed that COPD group samples have higher single-gene mRNA expression levels, with RNASE2 mRNA levels significantly higher than those in the normal control group (p=0.032). However, CD63 and ERAP2 did not show a difference between the two groups. Single-gene GSEA enrichment analysis showed that the 3 highly expressed differential genes in the COPD-smoking group are enriched in pathways related to allogeneic graft rejection, cobalamin transport and metabolism, and virus-related immune activation (Figure 4D–4F). Subsequent diagnostic ROC analysis revealed that the above 3 feature genes have high efficacy in diagnosing COPD (AUC: 0.586–0.727) (Figure 4G–4I). We then performed protein-protein interaction (PPI) analysis and using the STRING database to investigate the molecular interaction patterns of these genes. The results (Figure 4J and K) indicated that the primary interacting partners of CD63, RNASE2, and ERAP2 were CTSD, CTSA, and HLA-B, respectively. We also developed a clinical diagnostic array based on the three feature genes (Figure 4L).

Figure 4 Single-Gene Assessment, Protein-Protein Interaction NETwork Analysis, and Nomogram Construction in smoking COPD group. (A–C) present the gene expression group statistics for three smoking COPD NDEGs derived from machine learning on the GSE37768 dataset. The blue bars signify the normal control group, and the red bars represent the COPD group, with the y-axis indicating the level of gene expression. (D–F) depict the GSEA enrichment of differential genes associated with the high expression of the three single genes. (G–I) showcase the diagnostic dependency ROC curves for assessing the efficacy of single genes, with the yellow area corresponding to the area under the curve. (J and K) display the PPI protein-protein interaction NETworks and protein interaction diagrams based on STRING database analysis. (L) is a gene diagnostic Nomogram constructed from the three genes.

In the COPD-nonsmoking group, GRIPAP1, NHS, EGFLAM, and GLUL expression levels are all significantly elevated compared to the control group (P=0.0083, p=0.015, p=0.0024, and p=0.0016) (Figure 5A–D). The 4 highly expressed differential genes in the COPD-nons moking group were enriched in pathways related to IgA production, allogeneic graft rejection, apoptosis, and IL-17 signaling (Figure 5E–H). Subsequent diagnostic ROC analysis revealed that the above 4 feature genes have high efficacy in diagnosing COPD (AUC: 0.735–0.870) (Figure 5I-L). The protein interaction network for the four feature genes in the COPD-nonsmoking group comprises key players such as ENO1, GAPDH, and CD63, with comprehensive interaction maps presented in Figure 5M and N. A diagnostic array diagram for clinical application was developed based on the four feature genes (Figure 5O).

Figure 5 Single-Gene Assessment, Protein-Protein Interaction NETwork Analysis, and Nomogram Construction in nonsmoking COPD group. (A–D) present the gene expression group statistics for four nonsmoking COPD NDEGs derived from machine learning on the GSE37768 dataset. The blue bars signify the normal control group, and the red bars represent the COPD group, with the y-axis indicating the level of gene expression. (E–H) depict the GSEA enrichment of differential genes associated with the high expression of the three single genes. (I–L) showcase the diagnostic dependency ROC curves for assessing the efficacy of single genes, with the yellow area corresponding to the area under the curve. (M and N) display the PPI protein-protein interaction NETworks and protein interaction diagrams based on STRING database analysis. (O) is a gene diagnostic Nomogram constructed from the three genes.

Systematic Validation of Machine Learning Model Efficacy

Next, we conducted efficacy validation of the machine learning-constructed Models A and B across multiple datasets. Initially, analysis of the training set GSE37768 dataset showed that COPD group samples had significantly higher risk scores for both Model A and Model B (p<0.05); while the ROC curves indicated good diagnostic efficacy (AUC=71.388 and AUC=87.222), with Model B outperforming Model A, further distinguishing between smoking and nonsmoking samples (Figure 6A–D). Additionally, we used the external dataset GSE54837 for efficacy assessment. The results, as shown in Figure 6E–H, indicate that in Model A, both COPD stage high and stage low groups had higher scores than the smoking control group (p=0.044, p=0.043); while in Model B, only the COPD stage high score was significantly higher than the smoking control group (p=0.0097); Models A and B showed similar diagnostic efficacy (AUC=60.245 and AUC=62.287). Similarly, analysis of another independent COPD dataset GSE57148 revealed that, compared to the normal control group, COPD group samples had significantly higher risk scores in both Model A and Model B (p<0.001); while the ROC curves indicated that Model A had better diagnostic efficacy than Model B (AUC=70.228 and AUC=62.839) (Figure 6I–L).

Figure 6 Validation of the Risk Scoring Model with Internal and External Datasets. (A–D) depict the validation outcomes of the internal dataset GSE37768 using Model A and Model B. In (A and B), the y-axis denotes the computed risk scores, with violin plots in various colors representing different groups. (C and D) exhibit the diagnostic dependency ROC curves for Model A and Model B, respectively, with the surrounding hues indicating the confidence intervals. (E and F) and (I and J) illustrate the disparities in group model scores for the external validation datasets GSE54837 and GSE57148. (G and H) and (K and L) present the diagnostic dependency ROC curves for Model A and Model B on the external validation datasets GSE54837 and GSE57148, respectively, with the surrounding hues signifying the confidence intervals.

Distribution Patterns of RNASE2 and NHS and Protein Structure Prediction Based on Alphafold

Since the three machine learning methods assigned the highest importance coefficients to the RNASE2 and NHS genes, we conducted a systematic investigation and analysis of their distribution. The results, as shown in Figures Supplementary 7–8A, indicate that RNASE2 and NHS genes are expressed to varying degrees in multiple human organs, with RNASE2 being highly expressed in the spleen and lungs, and NHS in the retina and cervix; further analysis of their distribution in immune cells revealed that, in addition to high expression in neutrophils, RNASE2 and NHS genes are also highly expressed in PBMCs, monocytes, macrophages, basophils, and myeloid DC cells to varying degrees (Supplementary Figures 7, 8B and C). Subcellular localization revealed that RNASE2 exists in a membrane secretion manner, while NHS is involved in the construction of the cell membrane scaffold (Supplementary Figures 7–8D). Finally, we used Alphafold version 2.0 to construct protein structure models for RNASE2 and NHS (Supplementary 7–8E).

qRT-PCR and ELISA Validation of Alveolar Lavage Fluid Samples

Lastly, we extracted RNA from the sediment of alveolar lavage fluid samples from 10 COPD patients for RT-PCR validation and used the ELISA method to detect the relative concentration of MPO in the supernatant. The results, as shown in Figure 7A and B, indicate that the relative expression levels of RNASE2 and NHS genes in the alveolar lavage fluid of COPD stage-high patients are significantly higher than those in COPD stage-low patients. After constructing the MPO-ELISA standard curve, the results of converting sample OD values to relative concentrations using the formula y = 1220.5x - 118.1 showed that the supernatant of alveolar lavage fluid from COPD stage-high patients had a higher concentration of MPO (Figure 7C and D).

Figure 7 Quantitative RT-PCR and ELISA Analysis of RNASE2 and NHS Gene Expression and MPO Concentration in Clinical Samples. (A and B) display the outcomes of quantitative real-time polymerase chain reaction (qRT-PCR), where the red bars indicate COPD patients who need mechanical ventilation support (stage high), and the blue bars indicate those who do not (stage low). The y-axis denotes the relative expression levels of the target genes normalized against the GAPDH housekeeping gene, represented as the relative expression ratio (ΔCt value). (C) illustrates the ELISA standard curve that was established, while (D) shows the relative concentration of myeloperoxidase (MPO) calculated after calibration against the standard curve. The red bars in (D) represent COPD patients requiring mechanical ventilation support (stage high), and the blue bars represent those not requiring such support (stage low).

Discussion

This study innovatively combines single-cell sequencing with large-scale transcriptomics and employs three classic machine learning methods to identify NETs signature genes closely associated with poor prognosis in COPD. These genes were used to construct models for predicting outcomes in different COPD populations based on smoking differences. The effectiveness of these models was preliminarily validated using alveolar lavage fluid samples.

Inflammatory mechanisms play a decisive role in the development of COPD, characterized by increased and activated macrophages and neutrophils.¹⁴ An increase in sputum neutrophils is a significant feature of COPD, and this neutrophilic inflammation is induced by cigarette smoke, bacteria, viruses, and oxidative stress.¹⁰ Current research indicates that sputum from all levels of COPD patients, whether exacerbated or stable, is characterized by the presence of a large number of NETs and NET-forming neutrophils.^40,41 This is highly consistent with the results of our MCPcounter immune infiltration analysis study. Additionally, we detected varying concentrations of MPO in the alveolar lavage fluid of COPD patients. NETosis related to COPD is significantly associated with lung function impairment, which is considered one of the best indicators of disease severity, further strengthening the relevance and accuracy of our results. All characteristics of NETosis are more commonly observed in the sputum of exacerbated patients and in GOLD stage 3 and 4 patients.²⁰ Nevertheless, NETosis is clearly detected in the sputum of stable COPD and GOLD stage 1 and 2 patients.²⁰

The formation of NETosis can be understood as a result of neutrophil autophagy, an essential component of innate immunity.¹³ This programmed cell death through self-destruction is observed in numerous inflammatory diseases and cancers.^42,43 Their formation not only helps to confine pathogens and cancer cells but also aids in clearing pathogenic antigens due to the highly basic histones and degradative enzymes within their structure.⁴³ However, excessive neutrophil infiltration and activation leading to a large amount of pathological NETosis formation and cytokine storms may exacerbate lung function deterioration in patients with respiratory diseases.^17,44 Studies have shown that NETs formation is not limited to exacerbations but also exists in stable COPD and is correlated with the severity of airflow limitation.²⁰ NETs are major contributors to chronic inflammation and lung tissue damage in COPD most of the time. Pharmacologists have been working on manipulating NETs formation and administering DNase to degrade the NETs scaffold to reduce the pathogenic effects of NETs.⁴² With the development of genomics in recent years, a series of gene signatures related to NETs formation have been gradually identified in various diseases and have been well applied in predicting ischemia-reperfusion injury after lung transplantation, sepsis, and thrombotic diseases.^22,23 In our study, we used single-cell sequencing clustering analysis to identify differentially expressed genes in neutrophils and successfully identified NETosis-related genes associated with poor prognosis in COPD using machine learning combined with high-level gene signatures in COPD and ARDS from transcriptomics. Similarly, we detected higher levels of NETosis-related genes and MPO concentrations in alveolar lavage fluid samples from severe COPD patients requiring ventilator support.

A substantial body of research indicates that the genetic impact of smoking behavior increases over time.⁴⁵ Long-term smoking suppresses DNA repair mechanisms and causes genomic changes, leading to the occurrence and development of respiratory diseases and even cancer.⁴⁶ Non-smoking respiratory disease patients often have unique genetic biological characteristics and genetic susceptibility.⁴⁷ Interestingly, we found differences in gene transcriptomes between smoking and non-smoking COPD, leading us to construct separate NETosis-related gene models for predicting outcomes in smoking and non-smoking COPD patients. These models were preliminarily validated for effectiveness in smoking and non-smoking subgroups using multiple external datasets. Feature importance evaluation by the three machine learning methods showed that RNASE2 and NHS genes play a key role in smoking and non-smoking groups, respectively. Human RNase2 is a secreted protein expressed in white blood cells and is reported to have antiviral activity against single-stranded RNA viruses.⁴⁸ Current research also finds that RNase2 specifically cleaves cellular ncRNA.⁴⁹ During oxidative stress responses, Schlafen 2 (SLFN2) protein can protect tRNA from RNase cleavage.⁵⁰ The immunomodulatory role and potential targeting of RNase2 in human macrophage cell RNA populations have been preliminarily validated.⁴⁹ However, whether the high expression of the RNASE2 gene mediates neutrophil autophagy and further promotes NETosis formation during the inflammatory state of COPD exacerbations remains unconfirmed. NHS gene mutations are found to be associated with the X-linked hereditary disease Nance-Horan syndrome (NHS).^51,52 Early studies found extra expression of NHS in the kidneys, lungs, and thymus.⁵¹ Our subcellular localization found that NHS seems to be closely related to the scaffold structure of the cell membrane. There is currently very limited research on the NHS gene in immune cells.

In summary, the innovation of this study lies in the use of machine learning combined with single-cell omics and transcriptomics to screen for NET signature genes in smoking COPD and non-smoking COPD and to construct two gene models suitable for predicting outcomes in different COPD groups. However, the limitation of this study is the small number of clinical samples used for validation, which may cause some bias in the data. Additionally, there are no fixed indicators for the detection of neutrophil extracellular traps, often requiring multiple indicators and various detection techniques for identification. This study preliminarily uses MPO alone as a detection indicator for neutrophil extracellular traps, which can only serve as a preliminary reference, and the above shortcomings will be gradually improved in our subsequent research.

Conclusion

In conclusion, this study identified NET signature genes in smoking COPD and non-smoking COPD using various machine learning methods combined with single-cell omics and transcriptomics, and constructed COPD NDEGs Model A and Model B. The specificity of the feature genes RNASE2 and NHS and the good predictive efficacy of the scoring models were validated through internal training sets, external datasets, and clinical samples.