Pyroptosis-Related Genes as Diagnostic Markers in Chronic Obstructive Pulmonary Disease and Its Correlation with Immune Infiltration

Abstract

Background

Chronic obstructive pulmonary disease (COPD) stands as a predominant cause of global morbidity and mortality. This study aims to elucidate the relationship between pyroptosis-related genes (PRGs) and COPD diagnosis in the context of immune infiltration, ultimately proposing a PRG-based diagnostic model for predicting COPD outcomes.

Methods

Clinical data and PRGs of COPD patients were sourced from the GEO database. The “ConsensusClusterPlus” package was employed to generate molecular subtypes derived from PRGs that were identified through differential expression analysis and LASSO Cox analysis. A diagnostic signature including eight genes (CASP4, CASP5, ELANE, GPX4, NLRP1, GSDME, NOD1and IL18) was also constructed. Immune cell infiltration calculated by the ESTIMATE score, Stroma scores and Immune scores were also compared on the basis of pyroptosis-related molecular subtypes and the risk signature. We finally used qRT – PCR to detect the expression levels of eight genes in COPD patient and normal.

Results

The diagnostic model, anchored on eight PRGs, underwent validation with an independent experimental cohort. The area under the receiver operating characteristic (ROC) curves (AUC) for the diagnostic model showcased values of 0.809, 0.765, and 0.956 for the GSE76925, GSE8545, and GSE5058 datasets, respectively. Distinct expression patterns and clinical attributes of PRGs were observed between the comparative groups, with functional analysis underscoring a disparity in immune-related functions between them.

Conclusion

In this study, we developed a potential as diagnostic biomarkers for COPD and have a significant role in modulating the immune response. Such insights pave the way for novel diagnostic and therapeutic strategies for COPD.

Keywords:

• chronic obstructive pulmonary disease
• pyroptosis
• gene
• biomarkers
• immune Infiltration

Introduction

Chronic Obstructive Pulmonary Disease (COPD) is a significant medical concern worldwide due to its high incidence and mortality rates. Projections suggest that by 2030, COPD will rise as the third leading cause of death globally.¹ Research has shown that smoking is an important risk factor for patients with COPD, and harmful components in cigarette smoke can activate the innate immune system, recruit inflammatory cells to target lung tissue, and mediate the process of inflammation.² A retrospective study has found that smoking cessation can enhance the efficacy of single inhaler triple therapy, leading to better clinical and functional outcomes six months after commencement.³ Liu et al reported pyroptosis influences the development of COPD.⁴ Currently, there is evidence that cigarette smoke extract induces pyroptosis through the ROS/NLRP3/caspase-1 pathway.⁵ Therefore, it is particularly important to explore the regulatory mechanisms related to inflammation in COPD.

Pyroptosis is described proinflammatory programmed cell death and could play a significant role in a variety of biological systems.⁶ This process, driven by the GSDMD pyroptosis driver, forms non-selective pores in the cell membrane. This alteration subsequently induces cell swelling, rupture, and the eventual release of cellular contents.^7,8 GSDMD comprises of two domains, namely the N-terminal gasdermin-N and the C-terminal gasdermin-C, which are connected by a linking loop. The cleavage of GSDMD between these domains is adequate to induce pyroptosis. Being a pore-forming protein, GSDMD typically exists in a state of autoinhibition.⁹ Recently, upon acute exposure of mice to cigarette smoke (CS), NLRP3 facilitates the activation of GSDMD in bronchial epithelial cells and macrophages located in the bronchoalveolar space, presumably via the activation of the NLRP3 inflammasome.¹⁰ Notably, the initiation of pyroptosis releases inflammatory factors that can potentiate a robust inflammatory response, serving as precursors to certain diseases.¹¹ It has been reported that the nanotuner platform can cause pyroptotic cell death in various gasdermin-E-positive human cancers, leading to enhanced anti-tumour efficacy and minimized systemic side effects.¹² Given that inflammation is pivotal in the pathogenesis of COPD, the identification and elucidation of pyroptosis-related genes (PRGs) can furnish potential biomarkers for COPD diagnosis.

Our extensive literature review on pyroptosis culminated in the identification of 31 known PRGs.^13–15 These were subjected to meticulous scrutiny, with specific genes earmarked for subsequent evaluation. Our research integrated mRNA data from three distinct cohorts, synthesized the role of PRGs, and delved into their correlation with COPD-related immune cell infiltration. It was discerned that immune cell infiltration exhibited substantial variations across the two PRG modification patterns. These patterns exhibited a direct relation to immune checkpoints and overall immune activity. Building on this, we devised a pyroptosis scoring system tailored to pinpoint COPD patients based on PRG modification patterns. Our study is to identify and evaluate potential screening markers associated with COPD for early diagnosis and intervention of COPD. Our research was accepted by the medical ethics review of the Anqing Municipal Hospital (2022.No.61).

Materials and Methods

Source of Datasets

Datasets were sourced from the Gene Expression Omnibus (GEO) and annotated using R version 3.6.5,¹⁶ in conjunction with the GEO query package. These datasets detailed gene expression profiles and clinical data for patients diagnosed with COPD (Supplementary Table 1), all derived from Homo sapiens. COPD is heterogeneous, defining features include persistent airflow obstruction and respiratory symptoms. We primarily incorporated data from three distinct databases: (GEO dataset information summary Supplementary Tables 1–3).

GSE76925:¹⁷ Comprising 40 normal controls and 111 COPD patients (Supplementary Table 1).

GSE8545:¹⁸ Including 36 normal controls and 18 COPD patients (Supplementary Table 2).

GSE5058:¹⁹ Featuring 24 normal controls and 15 COPD patients (Supplementary Table 3).

For processing, annotation, and normalization of the datasets, we employed the hgu133plus2.db²⁰ and illuminaHumanv4.db²¹ annotations. As a result, the gene expression matrices for these datasets were duly acquired. In addition to these datasets, we conducted a comprehensive literature review pertinent to pyroptosis. Based on this review, we identified and shortlisted 31 recognized pyroptosis-related genes for subsequent analyses.^13–15

Express Mechanism of Pyroptosis-Related Genes

Firstly, we obtained an expression matrix of PRGs from the dataset GSE76925. Boxplot was used to visualize the expression of the PRGs between normal controls and COPD patients. Subsequently, the R programme RCircos²² was used to plot the circos plot for the differential expression across chromosomes. Expression pattern of PRGs was visualized with R function “pheatmap”.²³ Then, PPI networks of PRGs were constructed using STRING.²⁴ Correlations were visualized using the R corrplot package.²⁵ KEGG and GO term gene set enrichment analysis (GSEA) was done using clusterProfiler.²⁶

Modeling COPD Based on Pyroptosis-Related Genes

A logistic regression analysis was executed on PRGs of the dataset GSE76925 utilizing R, specifically leveraging the glm function.²⁷ Genes linked to the diagnosis of COPD were filtered based on variables with a significance level of p < 0.05. The LASSO²⁸ serves as a shrinkage estimation technique. This methodology involves the construction of a penalization procedure to compress numerous regression coefficients to zero, thus facilitating variable selection. Leveraging the pyroptosis genes associated with COPD diagnosis, a LASSO classifier was designed to pinpoint an optimized feature subset.

Risk scores for each sample across datasets GSE76925, GSE8545, and GSE5058 were computed based on the aforementioned model. ROC curves were subsequently formulated using the R package “pROC”.²⁹ A violin plot was employed to vividly depict the risk distribution between COPD and the control group across the datasets GSE76925, GSE8545, and GSE5058. To forecast the binding miRNA of COPD diagnostic markers, the miRDB³⁰ platform was utilized, culminating in the construction of mRNA-miRNA interaction networks. These networks were visualized employing the Cytoscape software.³¹

Construction of PRG Molecular Subclusters

Consensus clustering was performed using the ConsensusClusterPlus package in R to identify optimal profiles.³² Differentially expressed genes (DEGs) were identified using the limma package³³ in R. To select DEGs for in-depth analysis, we adopted the following criteria: a log2 fold change (FC) ≥1, and a false discovery rate (FDR) <0.05. The ESTIMATE R-package³⁴ was employed to compute both ESTIMATE, Stroma scores and Immune scores within the GSE76925 dataset. The ComplexHeatmap package³⁵ was utilized to generate heatmaps, illustrating the variation in clinical data distributions (such as Age, Sex, Race, BMI, FEV1, and FEV1/FVC subtype) across different groups. For molecular markers of immune-checkpoint subtypes, we selected CD274, CTLA4, HAVCR2, IDO1, LAG3, and PDCD1. Meanwhile, CD8A, CXCL10, CXCL9, GZMA, GZMB, IFNG, PRF1, TBX2, and TNF were chosen as molecular markers for immune-activity subtypes. Box plots were subsequently generated to visualize the immune-related characteristics of each subtype.

Gene Ontology /Kyoto Encyclopedia of Genes (GO/KEGG) Enrichment Analysis

The “GO plot” software package in R was used for GO³⁶ and KEGG³⁷ pathway enrichment analyses. GO included cell composition, biological processes, and molecular functions. Metascape was used to identify significantly enriched terms for GO and KEGG terms for differentially expressed genes (p value ≤0.01).

Gene Set Enrichment Analysis(GSEA)

Annotating potential pathway functional enrichments between gene sets using GSEA.³⁸ Gene set enrichment analysis of the gene expression matrix was used to analyze the gene expression profiles of different subtypes of GSE76925 using clusterProfiler. The “c2.cp.kegg.v7.4.symbols.gmt” was selected as a reference in the Molecular signatures database (MSigDB)³⁹ for GSEA. Statistical significance was defined as p≤0.1. The gseaplot application EnrichmentMap was used to visualize functional enrichment.

Weighted Correlation Network Analysis (WGCNA)

The network analysis was performed using an R package of WGCNA⁴⁰ in genes of GSE76925.Then, the soft threshold value of the correlation matrix was set to 4(R²=0.879) and height to 0.25 to choose parameter most relevant of COPD. With default parameters, Metascape was used to analyze Gene Ontology/KEGG terms and visualize enrichment results, identifying significant terms with a maximum overlap of 3 terms. p-value< 0.01, and a minimum enrichment of 1.5.

Immune Cell Infiltration Analysis and PRG Analysis

Gene expression profiles and immune cell files of GSE76925 were used to extract an immune cell infiltrate matrix in R using CIBERSORT,⁴¹ where the source code and the corresponding immune cell file were downloaded. The heatmaps for immune infiltration of 22 immune cells from distributions for each sample were plotted using the R package ”pheatmap”. Correlation heatmaps of immune cell infiltration were calculated in R and visualized with the corrplot package. We analyzed the correlation between pyroptosis related genes and immune cell infiltration and then visualized the results using the PHEATMAP package.

Human Clinical Specimens

RNA samples of COPD and normal person were obtained from Anqing Municipal Hospital, AnHui, P. R. China. The clinical information for these two groups was listed in Supplementary Table 4. The study protocol was approved by Anqing Municipal Hospital Research Ethics Committee. All patients signed clinical study consent forms. Written informed consent for publication of their details was obtained from the patient.

The diagnosis of COPD is based on the presence of risk factors, symptoms, and spirometry demonstrating persistent airflow obstruction. According to GOLD guidelines, post-bronchodilator testing with an FEV1 to forced vital capacity ratio less than 70% is recommended for the diagnosis of obstruction. People with COPD should receive smoking cessation, vaccinations, inhaler therapy, physical activity, pulmonary rehabilitation, long-term oxygen therapy and non-invasive ventilation.⁴²

Quantitative Reverse Transcription-Polymerase Chain Reaction (qRT-PCR)

Total RNA was extracted from whole blood samples of COPD patients and normal individuals using a combination of Trizol reagent (Invitrogen) and silica-based adsorption column method. Then the total RNA was reverse transcribed into cDNA using the HiScript 1st Strand cDNA Synthesis Kit (cDNA first-strand synthesis kit.).

Statistical Analysis

Analyses were performed using the R software (version 4.2.1). To evaluate the diagnostic impact of gene expression on normal and disease conditions, ROC curves were plotted and analyzed using the pROC package. Differential gene expression was assessed using the wilcox. test function in R. Significance levels were defined as: *p <0.05; **p <0.01; ***p <0.001; ****p<0.0001. Regarding the method of genetic identification, we performed univariate logistic regression analysis on the PRGs in GSE76925, selecting genes associated with COPD diagnosis based on a significance level of p<0.05.

Correlations between immune cells were quantified using the Pearson correlation coefficient test. The strength of these correlations was categorized as follows:

• Strong correlation: r≥ 0.8

• Medium correlation: 0.6≤ |r|≤ 0.79

• Moderate correlation: 0.4≤ |r|≤ 0.59

• Weak correlation: 0.2≤ |r|≤ 0.39

• No correlation: |r| < 0.2

Significant correlations between functional genes were visualized using the corrplot function from the corrplot package in R. Analyses examining associations between PRGs and immune cells were conducted using the psych package in R.

Results

Data Download and Pre-Processing

The mRNA expression profile datasets of GSE76925, GSE8545 and GSE5058 (Table 1) were downloaded from GEO (http://www.ncbi.nlm.nih.gov/geo/). Data processing of clinicopathological characteristics was performed with the GEOquery, Illuminahumanv4.db and HG U133p lus2.D B packages. Expression data and annotations of genes were then processed using MAS5.0 normalization with the “Affy” Bioconductor package. Gene annotation files were downloaded from HGNC. After matching the expression matrix contains 15,575 genes and 151 samples from GSE76925, 16,930 genes and 54 samples from GSE8545 and 16,930 genes and 39 samples from GSE5058. We searched the literatures related to pyroptosis and obtained 31 pyroptosis related genes (Table 2).

Table 1 Summary of Those 3 GEO Datasets of Chronic Obstructive Pulmonary Disease (COPD)

ID	GSE Number	Platform	Samples	Disease	References
1	GSE76925	GPL10558	111 patients and 40 controls	COPD	[17]
2	GSE8545	GPL570	18 patients and 36 controls	COPD	[18]
3	GSE5058	GPL570	15 patients and 24 controls	COPD	[19]

Abbreviation: COPD, chronic obstructive pulmonary disease.

Table 2 The 31 Original Pyroptosis-Related Genes That Were Used in This Study

Genes	Name
AIM2	Absent in melanoma 2
CASP1	Cysteine-aspartic acid protease-1
CASP3	Cysteine-aspartic acid protease-3
CASP4	Cysteine-aspartic acid protease-4
CASP5	Cysteine-aspartic acid protease-5
CASP6	Cysteine-aspartic acid protease-6
CASP8	Cysteine-aspartic acid protease-8
CASP9	Cysteine-aspartic acid protease-9
ELANE	Elastase, neutrophil expressed
GPX4	Glutathione peroxidase 4
GSDMB	Gasdermin B
GSDMC	Gasdermin C
GSDMD	Gasdermin D
GSDME	Gasdermin E
IL1B	Interleukin 1 beta
IL6	Interleukin 6
IL18	Interleukin 18
NLRC4	NLR family CARD domain containing 4
NLRP1	NLR family pyrin domain containing 1
NLRP2	NLR family pyrin domain containing 2
NLRP3	NLR family pyrin domain containing 3
NLRP7	NLR family pyrin domain containing 7
NOD1	Nucleotide binding oligomerization domain containing 1
NOD2	Nucleotide binding oligomerization domain containing 2
PJVK	Pejvakin/deafness, autosomal recessive 59
PLCG1	Phospholipase C gamma 1
PRKACA	Protein kinase cAMP-activated catalytic subunit alpha
PYCARD	PYD and CARD domain containing
SCAF11	SR-related CTD associated factor 11
TIRAP	TIR domain containing adaptor protein
TNF	Tumor necrosis factor

Expression Mechanism of Pyroptosis Related Genes

In this study, 31 PRGs derived from the literature were screened for analysis in GSE76925. Moreover, box plots show CASP4, CASP5, ELANE, GPX4, GSDMB, GSDME, NLRP1, NOD1 within highly expressed genes and CASP6, IL18 within lowly expressed ones as compared with the normal samples (Figure 1A). Chromosome map shows that CASP4, CASP5, IL18 locate in chr11; GSDME, NOD1 locate in chr7, and CASP6 locates in chr4 (Figure 1B). There is a distinct difference in the expression patterns of PRGs between COPD and normal samples through a heatmap (Figure 1C). An interactome network graph was generated using STRING (Figure 1D). Figure 1E depicts a heat map representation of correlation between pyroptosis related genes. GO analysis revealed that PRGs were associated with interleukin-1 beta production, interleukin-1 production, positive regulation of interleukin-1 beta production, inflammasome complex, cytosolic part, membrane raft and cysteine-type endopeptidase activity involved in apoptotic process, cysteine-type endopeptidase activity, cysteine-type peptidase activity, cysteine-type endopeptidase activator activity involved in apoptotic process (Figure 2A). KEGG pathway analysis showed PRGs were enriched in pathways, such as NOD-like receptor, TNF signaling pathway, IL-17 signaling pathway, Apoptosis and Necroptosis (Figure 2B). We used the R package Pathview to visualize pathway of apoptosis composing KEGG pathways (Figure 2C).

Figure 1 Expression mechanism of pyroptosis-related genes. (A) Boxplots of pyroptosis-related genes, Blue represents the normal group and red represents the COPD group (*p<0.05; **p<0.01; *** p<0.001; ****p<0.0001). (B) Chromosomal mapping of Differential pyroptosis genes. (C) Heatmap of pyroptosis-related genes, Blue represents the normal group, red represents the COPD group; (D) PPI analysis the differentially expressed pyroptosis-related genes, The larger the node, the darker the color, the more important the gene. (E) Correlation heatmaps of pyroptosis-related genes, The color represents the strength of the relatedness, with red indicating negative correlation, and blue indicating positive correlation.

Abbreviations: COPD, chronic obstructive pulmonary disease; PPI, protein–protein interactions.

Figure 2 Functional enrichment analysis of pyroptosis-related genes. (A) GO top enrichment function dotplot; (B) KEGG top enrichment function dotplot; (C) Apoptosis pathway diagram.

Abbreviations: GO, gene-ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Construction of PRG Molecular Subclusters

A total of 8 pyroptosis related genes (CASP4, CASP5, ELANE, GPX4, GSDME, IL18, NLRP1, NOD1) with significant diagnosis were selected using univariate logistic regression analysis combined with LASSO in GSE76925 (p<0.05) (Figure 3A–D). A diagnostic model based on eight diagnostic markers was constructed using multiple logistic regression. Subsequently, we applied the above model to diagnose for each sample in GSE76925, GSE8545 and GSE5058. ROC curve analysis showed that risk scores were a strong diagnostic ability with an AUC of 0.809 in GSE76925, 0.765 in GSE8545 and 0.956 in GSE5058 (Figure 4A–C). The violin diagram confirms that the model has the same diagnostic ability (Figure 4D–F). Next, miRDB was used to screen the targeted relationship between miRNAs and diagnostic markers of PRGs. In addition, 7 mRNA and 174 miRNA were obtained, diagnostic marker interaction network of miRNA-pyroptosis was constructed (Supplementary Figure 1).

Figure 3 Construction of a diagnostic model for COPD. (A) Univariate Logistic regression analysis. (B and C) LASSO logistic regression algorithm screening of COPD diagnostic markers. (D) Multivariate Logistic regression analysis.

Abbreviation: LASSO, logistic least absolute shrinkage and selection Operator.

Figure 4 Validation of COPD diagnostic model. (A) GSE76925 ROC analysis of risk score. (B) GSE8545 ROC analysis of risk score. (C) GSE5058 ROC analysis of risk score. (D) GSE76925 violin plots. (E) GSE8545 violin plots. (F) GSE5058 violin plots.

Molecular Typing Analysis Based on PRGs

Based on 31 genes related to pyroptosis, a consensus cluster analysis was performed using the ConsensusClusterPlus R package to classify GSE76925 COPD group into two gene clusters: gene cluster C1 and gene cluster C2 (Supplementary Figure 2).

The “limma” R package was used to screen 105 DEGs associated with pyroptosis, and the filter criteria were |logFoldchange (FC) | >1 and false discovery rate (FDR) <0.05, with 7 genes being upregulated and 98 genes downregulated in GSE76925 COPD group. Upregulated and downregulated genes are plotted by volcano plot (Figure 5A). The expression patterns and clinical characteristics of pyroptosis related genes were significantly different between the C1 and C2 samples through heatmap(Figure 5B and C).

Figure 5 Subtype analysis. (A) Volcano plot of DEGs in C1 and C2, Red represents upregulated gene, blue represents downregulated gene, gray represents indifference gene; (B) Heatmap of PRGs, Red represents C2 group, blue represents C1 group; (C) Heatmap of clinical characteristics, Red represents C2 group, blue represents C1 group; (D) Bee swarm plot of risk score, Red represents C2 group, blue represents C1 group; (E) Violin plot of immune score, Red represents C2 group, blue represents C1 group; (F) Boxplots of immune biomarker, Red represents C2 group, blue represents C1 group (*p<0.05; **p<0.01; ***p<0.001).

Abbreviations: DEGs, differentially expressed genes; PRG, pyroptosis-related genes.

The C1 group has higher risk scores than the C2 group by beeswarm plot in COPD diagnostic model scores (Figure 5D). We applied the Estimate algorithm from GSE76925 (COPD samples) to calculate various immune scores: ESTIMATE score, immune score and matrix score. The violin chart showed that each score of C1 group also exhibited significantly higher immune score than the C2 group (Figure 5E). The C1 group in immune-checkpoint (CD274, CTLA4, PDCD1) was significantly higher than the C2 group. And immune activity gene (CD8A, GZMA GZMB, PRF1) in the C1 group also showed significantly higher than C2 group, while TBX2 in the C1 group lower significantly than in the C2 group (Figure 5F).

GO/KEGG Enrichment Analysis

We performed gene enrichment analysis of 105 PRGs of differential genes using the Metascape online tool. The function was screened using p <0.01, minimum count of 3 and enrichment factor >1.5. Overall, the function of DEGs was significantly associated with inflammatory response, response to lipopolysaccharide, response to molecule of bacterial origin, cellular response to cytokine stimulus, leukocyte migration, myeloid leukocyte migration and cytokine receptor binding, receptor ligand activity, chemokine receptor binding, and signaling receptor activator activity. DNA-binding transcription activator activity, growth factor activity, transcription regulator complex, RNA polymerase II transcription regulator complex, extracellular matrix, external encapsulating structure, TNF signaling pathway, IL-17 signaling pathway, NF-kappa B signaling pathway, NOD-like receptor signaling pathway, MAPK signaling pathway (Figure 6) (Table 3). The detailed results were shown in Supplementary Table 5.

Table 3 GO/KEGG Enrichment Analysis

Category	GO	Description	Count	LogP
GO Biological Processes	GO:0006954	Inflammatory response	25	−20.85304906
GO Biological Processes	GO:0032496	Response to lipopolysaccharide	18	−16.26035766
GO Biological Processes	GO:0002237	Response to molecule of bacterial origin	18	−15.77215212
GO Biological Processes	GO:0009617	Response to bacterium	23	−15.17001703
GO Biological Processes	GO:0071345	Cellular response to cytokine stimulus	22	−14.37853773
GO Biological Processes	GO:0060326	Cell chemotaxis	14	−13.7640133
GO Biological Processes	GO:0050900	Leukocyte migration	14	−13.13173926
GO Biological Processes	GO:0071396	Cellular response to lipid	18	−12.81834099
GO Biological Processes	GO:0006935	Chemotaxis	18	−12.77465739
GO Biological Processes	GO:0042330	Taxis	18	−12.74569332
GO Molecular Functions	GO:0005126	Cytokine receptor binding	15	−13.41803475
GO Molecular Functions	GO:0005125	Cytokine activity	14	−13.00331778
GO Molecular Functions	GO:0048018	Receptor ligand activity	17	−11.91524193
GO Molecular Functions	GO:0030546	Signaling receptor activator activity	17	−11.8170326
GO Molecular Functions	GO:0030545	Signaling receptor regulator activity	17	−11.1824344
GO Molecular Functions	GO:0042379	Chemokine receptor binding	7	−8.255128539
GO Molecular Functions	GO:0008009	Chemokine activity	6	−7.722573377
GO Molecular Functions	GO:0001664	Gprotein-coupled receptor binding	9	−6.108989859
GO Molecular Functions	GO:0005539	Glycosaminoglycan binding	10	−7.975209334
GO Molecular Functions	GO:0008201	Heparin binding	7	−5.64462679
GO Cellular Components	GO:0005667	Transcription regulator complex	13	−7.728286783
GO Cellular Components	GO:0090575	RNA polymerase II transcription regulator complex	6	−3.684595276
GO Cellular Components	GO:0031012	Extracellular matrix	7	−2.414412546
GO Cellular Components	GO:0030312	External encapsulating structure	7	−2.410283073
GO Cellular Components	GO:0044853	Plasma membrane raft	3	−2.139577205
KEGG Pathway	hsa04668	TNF signaling pathway	14	−17.55473142
KEGG Pathway	hsa04657	IL-17 signaling pathway	13	−16.92025671
KEGG Pathway	hsa04060	Cytokine-cytokine receptor interaction	16	−14.18269091
KEGG Pathway	hsa04064	NF-kappa B signaling pathway	9	−9.960989991
KEGG Pathway	hsa04621	NOD-like receptor signaling pathway	10	−9.015528688
KEGG Pathway	hsa05323	Rheumatoid arthritis	8	−8.886205988
KEGG Pathway	hsa04061	Viral protein interaction with cytokine and cytokine receptor	8	−8.633101909
KEGG Pathway	hsa05417	Lipid and atherosclerosis	9	−7.179598211
KEGG Pathway	hsa05146	Amoebiasis	7	−7.150791148
KEGG Pathway	hsa05144	Malaria	8	−11.11086597

Abbreviations: GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Figure 6 GO/KEGG enrichment analysis. (A) Network diagram of the top 20 enrichment functions the top 20 enriched terms in KEGG analysis. Cluster ID to indicate the color. Each node is an enriched term. (B) Network diagram of the top 20 enrichment functions. P value indicates the color. Each node is an enriched term. (C) GO BP enrichment results dot plot. The size of the dot represents the count number, and the color represents P value. (D) GO MF enrichment results dotplot. (E) GO CC enrichment results dot plot; (F) KEGG enrichment results.

Abbreviations: GO, gene-ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; BP, biological processes; MF, Molecular Function; CC, Cellular Components.

GSEA

We conducted GSEA to analyze the pathway enriched in expression profile of GSE76925 /COPD according to C1 and C2 using clusterProfiler. A p-value<0.05 was used as the threshold to determine the significant enrichment of KEGG pathways. (Supplementary Figure 3, Table 4). GSEA pathway of COPD in GSE76925 enrichment results showed that Salmonella infection, MAPK signaling pathway, Cytokine-cytokine receptor interaction, Human cytomegalovirus infection, NF-kappa B signaling pathway, PI3K-Akt signaling pathway, HIF-1 signaling pathway, Chagas disease, Pathways in cancer, Transcriptional misregulation in cancer, Apoptosis, Th17 cell differentiation, AGE-RAGE signaling pathway in diabetic complications, Coronavirus disease COVID-19, MicroRNAs in cancer, TNF signaling pathway, JAK-STAT signaling pathway markedly enriched (Supplementary Table 6) (Table 4).

Table 4 GSEA

Category	Enrich Function	NES	EnrichmentScore	Pvalue
KEGG Pathway	Salmonella infection	−1.712018306	−0.409682221	0.001876173
KEGG Pathway	MAPK signaling pathway	−1.841890851	−0.437790808	0.001886792
KEGG Pathway	Cytokine-cytokine receptor interaction	−2.594083879	−0.625896563	0.001904762
KEGG Pathway	Human cytomegalovirus infection	−2.120409696	−0.518117418	0.001915709
KEGG Pathway	NF-kappa B signaling pathway	−2.500995015	−0.667940292	0.001926782
KEGG Pathway	PI3K-Akt signaling pathway	−1.828888664	−0.427287909	0.001930502
KEGG Pathway	HIF-1 signaling pathway	−1.818057434	−0.485601113	0.001934236
KEGG Pathway	Chagas disease	−2.296923224	−0.617713744	0.001934236
KEGG Pathway	Pathways in cancer	−1.694577382	−0.3761112	0.001934236
KEGG Pathway	Transcriptional misregulation in cancer	−1.858152784	−0.464047357	0.001934236
KEGG Pathway	Prostate cancer	−1.846604718	−0.49660916	0.001934236
KEGG Pathway	C-type lectin receptor signaling pathway	−2.382461064	−0.638791792	0.001941748
KEGG Pathway	Tuberculosis	−2.228464932	−0.558434209	0.001941748
KEGG Pathway	Th17 cell differentiation	−1.909515919	−0.512727104	0.001945525
KEGG Pathway	T cell receptor signaling pathway	−2.045880988	−0.549342702	0.001945525
KEGG Pathway	AGE-RAGE signaling pathway in diabetic complications	−2.245180649	−0.602856965	0.001945525
KEGG Pathway	Coronavirus disease - COVID-19	−1.894609174	−0.462941453	0.001949318
KEGG Pathway	Parathyroid hormone synthesis, secretion and action	−1.790883333	−0.482644158	0.001953125
KEGG Pathway	Cellular senescence	−2.112588758	−0.535582749	0.001960784
KEGG Pathway	Human immunodeficiency virus 1 infection	−1.610105286	−0.398110802	0.001960784

Abbreviations: GSEA, Gene Set Enrichment Analysis;KEGG, Kyoto Encyclopedia of Genes and Genomes.

WGCNA

We performed robust weighted gene co-expression network analysis to expression profile in GSE76925/COPD. With the soft threshold set at 4 (scale-free R2 = 0.879) and the cut height at 0.25, we finally identified sixteen modules relevant to this subtype of COPD (Figure 7A–D). The association between modules and clinical characteristics was measured by the correlation between modular feature gene values and clinical features. Figure 7C and D showed that black module was most significantly associated with subtype of COPD. (correlation = 0.72, p = 4.6E-70). We uploaded black module into STRING database to build the protein–protein interaction (PPI) network. The sub-network in the top ranked-score were identified with the MCODE to obtained a module (Figure 7E). We considered that this module might play a special role in the diagnostic mechanism of COPD subtypes, then put this module gene into Metascape. The function was screened by p<0.01, minimum count of 3 and enrichment factor >1.5. The results showed differentially expressed genes are primarily related to NF-kappa B signaling pathway, Epstein-Barr virus infection, Human T-cell leukemia virus 1 infection, MAPK signaling pathway, TNF signaling pathway, IL-17 signaling pathway, inflammatory response, Coronavirus disease-COVID-19, cellular response to biotic stimulus, Chemokine signaling pathway (Figure 7F).

Figure 7 Network analysis identifies modules of co-expressed genes across disease. (A) Sample dendrogram and trait heat map; (B) Analys is of the scale-free fit index (left) and the mean connectivity (right) for various soft-thresholding powers; (C) The correlation of heatmap between subtype and module signature genes in COPD. (D) The correlation of scatterplot between subtype and black module. (E) The sub-networks with the highest of MCODE scores were identified, that is an important functional sub-networks in the black module. (F) GO\KEGG enrichment analysis of functional sub-networks.

Abbreviations: COPD, chronic obstructive pulmonary disease; GO, Gene Ontology.

Analysis of Correlation Between Immune Cell Infiltration and PRGs

Correlation analysis is results show (Figure 8A) that T cells follicular helper was positively correlated with GSDMB, GSDMD, NOD1, PLCG1, negatively correlated with CASP6. Macrophages M2 was significant negative correlated with GPX4, and positively with IL18. There was significant positively correlation between T cells CD4 memory resting and CASP6, Plasma cells and NLRP7, AIM2, T cells gamma delta was significant negatively correlated with IL18. Heat maps of 22 immune cells shows a clear negative correlation between Mast cells resting and Mast cells activated; between T cells CD4 memory resting and T cells follicular helper (Figure 8B). Heatmap of immune cell infiltration showed significant differences between COPD and the normal group (Figure 8C). Compared to normal group, COPD patients had higher proportions of immune cells with CD4 memory activated and CD4 memory resting, and there was a significant difference between COPD and the normal group (Figure 8D). Finally, the box plot of immune cell infiltration presents T cells CD4 memory resting and Macrophages M2 showed a higher level of immune infiltration compared to normal, whereas Plasma cells, T cells follicular helper and T cells gamma delta with a lower level (Figure 8E).

Figure 8 Correlation analysis between immune cell infiltration and pyroptosis. (A) Correlation analysis between pyroptosis genes and 22 immune cell infiltration, red represent positive correlation, blue represent negative correlation; (B) The correlation heat map of 22 immune cells infiltration, red represent negative correlation, blue represent positively correlation, The darker the colour, the stronger the correlation; (C) Heatmaps of immune cells infiltration between COPD group and normal group; (D) Box plots of immune cell infiltrate composition, (E) Box plots of two immune cell infiltrate composition, red represents COPD group; blue represents normal group, *p<0.05; **p<0.01; ***p<0.001.

Expression of PRGs in COPD and Normal Control

Quantitative real-time PCR (qRT-PCR) was performed to detect PRGs expression. The results showed that NOD1 and GSDME expression in COPD were significantly higher than in healthy individual (Figure 9A and B). However, expression of the six PRGs was downregulated in COPD compared with normal individual (Figure 9C–H).

Figure 9 Continued.

Figure 9 RT-qPCR analyze the mRNA of PRGs in blood between normal and COPD. (A) NOD1 mRNA (B) GSDME mRNA (C) IL-18 mRNA (D) CASP4 mRNA (E) CASP5 mRNA (F) ELANE mRNA (G) GPX4 mRNA (H) NLRP1 mRNA. *p<0.05; ***p<0.001; RT-qPCR: Real-time quantitative PCR.

Discussion

By 2030, COPD is projected to rank as the third primary cause of morbidity and mortality worldwide.¹ Current therapeutic strategies for COPD emphasize minimizing disease activity, mitigating symptoms, stalling disease progression, and ultimately reducing disease severity. A significant challenge in the management of COPD remains the accurate diagnosis and prognosis.

Recently, a novel method of programmed cell death, termed Pyroptosis, has come into the spotlight. Investigations have linked pyroptosis with tumor immunity and its therapeutic implications. Chen et al⁸ have showcased that nanoparticles can induce pyroptosis predominantly through a mitochondrial apoptosis pathway, triggered by PLC in early endosomes. Meanwhile, Chai et al⁴³ revealed the pivotal role of GSDMD-induced inflammatory cytokine release and pyroptosis in fostering host immune responses against infections. Meng et al¹³ provided evidence suggesting that pyroptosis scores could predict the prognosis and immunotherapy response in patients with cutaneous melanoma. To our knowledge, this study is the pioneer in exploring the role of Pyroptosis-related genes (PRGs) in COPD diagnosis.

Previous research has identified diagnostic biomarker signatures. Meng et al¹³ pinpointed an optimal panel of PRGs that could guide immunotherapy candidates and therapeutic strategies, anchored in comprehensive evaluations of tumor PRG modifications. By employing pyroptotic patterns, one can ascertain patient immunophenotypes and tailor immunotherapy approaches.¹⁴ Furthermore, Chao et al⁴⁴ introduced a prognostic model associated with PRGs in glioma patients. This model’s AUC values in the test set reached 0.669, 0.713, and 0.709 for 1-year, 3-year, and 5-year intervals, respectively. In contrast, our study identified eight PRGs with an AUC ranging from 0.765 to 0.956. Our findings suggest that a select combination of biomarkers is instrumental in diagnosing COPD.

Notably, CD8+ T cells and NK cells have been evidenced to instigate tumor pyroptosis via the GSDMB-GZMA axis, driven by interferon-γ (IFN γ). These immune cells could potentially target GSDMB-expressing cancer cells to enhance antitumor immunity.⁴⁵ Additionally, the nuclear PD-L1 family, when overexpressed, could shift TNF-induced apoptosis to pyroptosis, leading to tumor necrosis and potentially promoting tumor growth.⁴⁶ These studies underscore the intricate ties between pyroptosis and immune responses in tumor treatments. Yet, the domain of COPD remains untouched in the context of PRG bioinformatics.

Our research endeavors led us to devise a pyroptosis score system, gauging the diagnostic potential of PRGs in individual COPD patients. Based on the expression of 8 PRGs, we formulated two distinct pyroptosis modification paradigms. Our diagnostic classifier integrating these PRGs showcased robust discrimination in the training sets (AUC values: GSE76925=0.809, GSE8545=0.765) and held up impressively in the test set GSE5058 (AUC=95.6%). Beyond identifying the eight key PRGs (CASP4, CASP5, ELANE, GPX4, GSDME, IL18, NLRP1, NOD1) as potential biomarkers for COPD, we have also sculpted a diagnostic model rooted in these PRGs. Both internal and external validations reinforced the model’s efficacy in distinguishing COPD patients from control groups.

In addition, our study found five PRGs (NOD1, GSDME, GPX4, NLRP4, CASP4) were significantly different between COPD and normal control, which was consistent with the results of the box plot analysis of PRGs in GSE76925. NOD1 and GSDME were overexpressed in the COPD, while GPX4, NLRP1, CASP4 is underexpressed in the COPD. This finding is inconsistent with the results reported in the literature, which may be due to differences in the clinical samples collected.

The cytoplasmic nucleotide-binding oligomerization domain-1 (NOD1), a member of the NLR family, is widely expressed in various cell types.^47,48 The N-terminal of NOD1 contains caspase activation/recruitment domains (CARDs), which are highly expressed in inflammatory cells. The activation of NF- κ B signaling pathway leads to the secretion of inflammatory cytokines.^49,50 In advanced COPD stages, the upregulation of TLR4 and NOD1 in bronchial epithelia is connected with heightened inflammation and an increased bacterial load of P. aeruginosa.⁵¹ Elevating NOD1 expression in COPD mouse models exacerbated lung function impairments caused by cigarette smoke, as manifested by increased airway resistance, diminished functional residual capacity, and more pronounced pulmonary injury.⁵² In our study, T cells follicular helper displayed a positive correlation with GSDMB, GSDMD, and NOD1.

GSDME belongs to the Gasdermin family, which also includes GSDMA, GSDMB, GSDMC, GSDMD, and DFNB59.⁹ Gasdermin E possesses a autoinhibitory dual-domain architecture, comprising the N-terminal and C-terminal domains of Gasdermin.⁵³ As vital substrates for pyroptotic cell death, GSDME and GSDMD facilitate distinct forms of pyroptosis. GSDME-dependent pyroptosis is an inflammasome-independent pyroptosis without the assembly of inflammasomes. However, GSDME-N cleaved by caspase-3 has the pore-forming property to promote the release of cytokines such as IL-1β and IL-18, which have been confirmed in most pathological and pharmacological studies.^54–56 Recently, using GSDME-deficient mice and human tubular epithelial cells (TECs), Xia et al found that cisplatin treatment led to pyroptosis in TECs, upregulation of GSDME-N expression, and induction of cell pyroptosis through pharmacological and genetic interventions. Silencing GSDME attenuated acute kidney injury and inflammation in mice. These findings suggest that Caspase-3/GSDME-induced pyroptosis and inflammation contribute to the development of acute kidney injury (AKI).⁵⁶ Gan et al demonstrate that GSDME plays a role in the development of Crohn’s disease (CD) by inducing GSDME-mediated pyroptosis, which leads to the release of proinflammatory cytokines.⁵⁷ Furthermore, knockout of GSDME switched the manner of death of A549 and human primary alveolar epithelial cells from pyroptosis to apoptosis upon H7N9 virus infection, and GSDME knockout mice survived H7N9 virus lethal infection.⁵⁸ GSDME also plays a significant role in Corona Virus Disease 2019 (COVID-19). Specifically, the activation of the caspase-3/GSDME axis triggers the demise of airway epithelial cell during infection with SARS-CoV-2, particularly when GSDMD is absent.⁵⁹ Both our study and the box plot analysis of PRGs in relevant literature consistently revealed a significant overexpression of GSDME in the COPD group compared to the control group.

GPX4 serves as a crucial regulator of ferroptosis, and it is capable of inducing ferroptosis in mouse tumor xenografts.⁶⁰ Recent findings indicate that GPX4 exhibits dual functionality in OA, serving both to regulate ferroptosis or oxidative stress and to modulate ECM degradation via the MAPK/NFkB signaling pathway.⁶⁰ Copper directly binds to GPX4 protein, leading to the formation of GPX4 aggregates and subsequent autophagic degradation of GPX4.⁶¹ Our research corroborated that M2 Macrophages were more prevalent in the COPD group, exhibiting a strong inverse correlation with GPX4 and a positive one with IL-18. In contrast, T cells gamma delta had a significant negative correlation with IL-18.

NLRP1, a member of the NOD-like receptor (NLR) family, possesses a domain known as FIIND (function-to-find domain).⁶² The NLRP1 inflammasome serves as an innate immune sensor for the SARS-CoV-2 3CL protease.⁵⁹ Inflammasome-dependent caspase-1 activity can result in pyroptosis.⁶³ CASP4/5 bind intracellular lipopolysaccharide (LPS), the primary structural component of the outer membrane of Gram-negative bacteria, with high specificity and affinity via their N-terminal CARD domain, triggering caspase activation and resulting in pyropoptotic cell death.⁶⁴ Inhibiting CASP4 results in compromised cell migration and adhesion to the extracellular matrix in epithelial cancer cells.⁶⁵ Alveolar macrophages isolated from asthma patients display elevated expression of caspase-4, which is the human homologue of caspase-11. This upregulation of caspase-4/11 suggests a significant role in allergic airway inflammation.⁶⁶

Additionally, our thorough examination revealed significant differences in immune cell infiltration patterns in relation to various pyroptosis types and their clinical implications. Notably, the GSEA result was significantly enriched in a myriad of pathways, including the MAPK signaling pathway, NF-kappa B signaling pathway, and HIF-1 signaling pathway, among others.

We employed CIBERSORT for an exhaustive analysis of COPD immune infiltration to delve deeper into the role of immune cell infiltration in this disease. We found the proportion of CD4 memory-resting T cells and M2 macrophages in COPD patients was significantly elevated compared to controls, while the presence of Plasma cells, T cells follicular helper, and T cells gamma delta was considerably lower.

Research indicates that pyroptosis affects immune system components, such as macrophages and dendritic cells (DCs).⁶⁷ Pyroptosis mediates interactions between innate and adaptive immunity, shaping the cancer microenvironment to induce an immunostimulatory response.⁶⁸ In our study, we explored the relationship between several pyroptotic genes and immune cell infiltration. We found these PRGs to be significantly correlated with immune cells, suggesting an interplay between pyroptosis and the immune response in COPD. The underlying molecular mechanisms remain to be elucidated.

Our study has certain limitations. Although there were differences in PRGs between the two groups in our clinical specimens and database data analysis, CASP4, CASP5, ELANE, GPX4, and NLRP1 were found to be expressed at lower levels in the COPD group compared to the control group in our clinical whole blood specimens. This finding is inconsistent with the results reported in the literature, which may be due to differences in the clinical samples collected. Our external validation was performed on two datasets with a relatively small sample size in the external validation set. Although the model’s AUC demonstrated reasonable discriminatory power, there is room for improvement. Consequently, our goal in this study was to identify and evaluate potential screening biomarkers associated with COPD for early diagnosis and intervention. Assessing their diagnostic accuracy and sensitivity in a larger population is a crucial direction for future research.

Conclusion

In this study, we identified eight potential PRGs—namely, CASP4, CASP5, ELANE, GPX4, GSDME, IL18, NLRP1, and NOD1—as potential diagnostic biomarkers for COPD. Notably, This is the first time these PRGs have been associated with COPD diagnosis. We conducted both internal and external validations of a PRGs diagnostic model for COPD using LASSO logistic regression. The AUC value of this model in the test set yielded values of 0.765, 0.809, and 0.956 for datasets GSE8545, GSE76925, and GSE5058, respectively.

To further investigate the role of the immune response in COPD, we employed CIBERSORT to assess immune cell infiltration using our gene signatures. Preliminary correlation analysis suggests that pyroptosis may be implicated in the immune response observed in COPD patients. Thus, our gene signature serves as potential markers for the screening of COPD. The insights gained from this study not only enrich our understanding of the pyroptosis genes in patients with COPD but also aim to screen for biomarkers to achieve early intervention, reduce the disease burden, and improve the quality of life for patients with COPD.

Ethical/Copyright Corrections

Medical ethics review (2022) No. 61 by Anqing municipal hospital.

Consent for Publication

We give our consent for publication.

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.

Disclosure

The authors report no conflicts of interest in this work.

Acknowledgment

We thank (https://www.helixlife.cn/) for its linguistic assistance during the preparation of the manuscript. We would like to acknowledge Yingchun Ji’s contribution to data collection for this article.

Data Sharing Statement

The data that support the findings of this study are openly available in GEO. We used three databases: GSE76925 (doi:10.1038/srep44232), GSE8545 (doi:10.1378/chest.07-2245), and GSE5058 (doi:10.1164/rccm.200705-795OC).

Additional information

Funding

This study was supported by the Health Research Program of Anhui (AHWJ2022b099 to S.H.M.) and Anhui Provincial Key Research and Development Project (2022e07020052 to Z.H.P).