We retrieved the gene expression data for BPD and control samples from the Gene Expression Omnibus (GEO) database. GSE32472 provides microarray profiles of blood samples from neonates with BPD, including microarray analyses of gene expression on the 5th, 14th, and 28th days of life. On day 28, after a more definitive diagnosis of BPD was established, 100 blood samples were chosen at random to verify the reliability of the selection process. These samples were collected from 38 controls and 62 patients with BPD. We calculated the Effect size d to evaluate the statistical power of the sample size via G*Power software (Faul et al., 2009), which confirmed the sample size was reasonable. The gene set enrichment analysis (GSEA) public database was searched with the term “pyroptosis.” Two gene sets were obtained, GOBP PYROPTOSIS and REACTOME PYROPTOSIS. In total, 40 PRGs were discovered. In addition, 33 PRGs were identified based on previous reports. Finally, 57 PRGs were identified (Supplementary Table S1).

Expression data for 55 PRGs were initially retrieved from the GSE32472 dataset for control and BPD samples. Afterward, we run the Wilcoxon test in the R software to identify the differentially expressed PRGs (DE-PRGs) between the two samples. A p-value less than 0.05 was chosen to serve as the criterion for statistical significance.

The Metascape database (http://metascape.org) was employed to explore the biological significance of the differentially expressed genes (Zhou et al., 2019). Gene Ontology (GO) and Kyoto Encyclopedia of Genomes (KEGG) pathway enrichment analyses were performed using specific genes. The minimum overlap required for statistical significance was set at three, and the p-value was less than 0.01.

With the aid of the glmnet package, the LASSO algorithm was applied to the data to minimize its dimensionality. To select markers and identify gene biomarkers for BPD, the DE-PRGs between control and BPD samples were retained. The SVM package in R was used to develop support vector machine-recursive feature elimination (SVM-RFE) models, and the average misjudgment rates of their 10-fold cross-validation were compared. In addition, optimum gene biomarkers for BPD were discovered using biomarkers produced from both algorithms that overlapped. The diagnostic capabilities of the optimal gene biomarkers were evaluated by generating the receiver operating characteristic (ROC) curve and determining the area under the curve (AUC), sensitivity, specificity, and accuracy. Also, the prediction function in the R package was used for the purpose of developing a logistic regression model based on the five marker genes. This was carried out to predict the types of samples present in the GSE32472 dataset. We used ROC curves to assess the logistic regression model’s diagnostic power.

The GSEA R package was used to perform this analysis. To further investigate the linked pathways of the five marker genes, the correlation between these marker genes and the other genes in the GSE32472 dataset was evaluated. The genes were then ranked, from highest to lowest, based on their correlations, and the ranked genes were referred to as the test set. The KEGG pathway set served as a predetermined set to discover its enrichment in the gene set.

The GSVA R package was used for this analysis. The KEGG pathway enrichment set was used as the background gene set to analyze each marker gene. Using the limma program, we simultaneously performed an analysis of the difference in GSVA scores that existed between groups that had high and low-expression levels of the marker genes. The screening condition difference were p < 0.05 and |t| > 2. When t was more than zero, the pathway was considered to be enriched in the group that had high expression levels; conversely, when t was less than zero, the pathway was considered enriched in the group that had low expression levels.

CIBERSORT is a technique for determining the cell fraction of complex tissues based on gene expression profiles (Newman et al., 2015). We used CIBERSORT to predict the proportion of 22 different types of infiltrating immune cells in each sample from the GSE32472 dataset. All immune cell type fractions were summed up to one for each sample.

mRNA–miRNA interaction pairs based on the five marker genes were predicted using the starBase database (Li et al., 2014). RNA sequences for the five marker genes were retrieved from the NCBI website, while human miRNA sequences were retrieved from the miRBase database. mRNA–miRNA nucleic acid binding was predicted using miRanda, and the binding score cutoff value was elevated from 140 to 170. The predicted miRNA was then searched in starBase and miRNA-lncRNAs were screened to obtain the ceRNA network of mRNA-miRNA-lncRNA.

A diagnostic nomogram based on selected candidate PRGs was established using the “rms” package in the R software. The consistency between our predicted values and reality was assessed using a calibration curve. Clinical impact curve analysis and decision curve analysis (DCA) were carried out to determine whether or not the decisions made based on the model were favorable to the patients.

The JCRB Cell Bank provided the human lung adenocarcinoma epithelial cell line, A549. It was grown in DMEM with 1% penicillin/streptomycin solution and 10% fetal bovine serum (Thermo Fisher Science, Waltham, MA). Experiments on hyperoxia were carried out in a Modular Incubator Chamber (MIC-101, Billirups-Rothenberg Inc., Del-Mar, CA) in accordance with the manufacturer’s instructions. In details, the hyperoxia group was cultivated in a MIC-101 incubator with 85% O2 and 5% CO2, whereas the normoxia (NO) group continued to be cultured in an incubator with 21% O2 and 5% CO2. The RNeasy Mini Kit was used to extract total RNA (Qiagen, Hilden, Germany). Utilizing ReverTra Ace, RNA was reverse transcribed into cDNA (Toyobo, Osaka, Japan). Using the LightCycler 480 (Roche Applied Science, Indianapolis, IN), quantitative real-time PCR was carried out using SYBR Premix Ex Taq II (Takara Bio, Shiga, Japan). Primers used were designed in PrimerBank (http://pga.mgh.harvard.edu/primerbank/) and synthesized by Generay (Shanghai, China). The sequences were as follows: NLRC4, forward 5′-TGC​ATC​ATT​GAA​GGG​GAA​TCT​G-3′ and reverse 5′- GAT​TGT​GCC​AGG​TAT​ATC​CAG​G-3′; CHMP7, forward 5′- AAG​CCT​CTC​AAG​TGG​ACT​CTT-3′ and reverse 5′- ACA​GAC​GAT​ACA​CCT​CCT​CAG-3′; NLRP2, forward 5′- TGG​CCT​GGA​GAT​AGC​AAA​GAG -3′ and reverse 5′- CAC​CAC​CGT​GTA​TGA​GAA​GGG -3′; NLRP6, forward 5′- CCT​ACC​AGT​TCA​TCG​ACC​AGA-3′ and reverse 5′- CTC​AGC​AGT​CCG​AAG​AGG​AA -3′; NLRP9, forward 5′- TTT​GGC​TTG​TTG​TGG​TAT​CTG​AA -3′ and reverse 5′- CTG​GGT​AAT​GTT​TGT​CCA​GCA-3′; and GAPDH, forward 5′- GGA​GCG​AGA​TCC​CTC​CAA​AAT-3′ and reverse 5′- GGC​TGT​TGT​CAT​ACT​TCT​CAT​GG -3′.

The BPD vs. control groups were compared by the Wilcoxon test. The relationship between the DE-PRGs was determined using Pearson’s correlation analysis. The Jvenn package was employed to generate the Venn diagrams. The ceRNA network was visualized using Cytoscape version 3.9.1. A p-value less than 0.05 was considered significant. The R package was used for all of the analyses that were conducted.

In total, 25 out of 55 PRGs, including nine downregulated and 16 upregulated genes, were found to be differentially expressed between control and samples in the GSE32472 dataset (Supplementary Table S2). The expression patterns of the DE-PRGs were visualized in a clustering heatmap (Figure 1A). Figure 1B presents the correlation between these genes. As shown in Figure 1C, D, NLRP6, NLRC4, NLRP9 and other 13 DE-PRGs were upregulated in BPD group, while CHMP7, NLRP2, and other 4 DE-PRGs were downregulated, further confirming the results in Figure 1A.

GO and KEGG pathway enrichment analyses were performed using Metascape. The top GO and KEGG pathways were pyroptosis, regulation of proteolysis, shigellosis, positive immune response regulation, legionellosis, and regulation of interleukin-18 production, indicating that DE-PRGs may be involved in the pathogenesis of BPD (Figure 2).

We then determined the diagnostic potential of the DE-PRGs. We then used the SVM-RFE and LASSO algorithms to screen the DE-PRGs in the GSE32472 dataset to differentiate BPD from control samples. Six BPD-related markers were selected using the LASSO regression with penalty parameter tuning carried out by 10-fold cross-validation (Figures 3A, B). The 25 DE-PRGs were then filtered using the SVM-RFE algorithm to determine the optimal combination of feature genes. The optimal feature genes were ultimately found to be 18 in number (maximal accuracy = 0.75, minimal RMSE = 0.25) (Figures 3C, D). Five marker genes (NLRP6, CHMP7, NLRP2, NLRC4, and NLRP9) were identified for further analysis after the marker genes were subjected to an intersection (Figure 3E). To further confirm the roles of the five marker genes in BPD, A549 cells were cultured in a special incubator containing 85% O2 and 5% CO2 to simulate the BPD state, while control cells were cultured in 21% O2 and 5% CO2. The RT-qPCR results in Supplementary Figure S1 showed that NLRP6, NLRC4, and NLRP9 were elevated in the hyperoxia group, while CHMP7 and NLRP2 were downregulated. The difference in the expression levels of the five genes between hyperoxia and normoxia groups showed the same trends as those between BPD and control samples.

A regression model based on the five marker genes was established using the glmnet package in R. Following this, ROC curves revealed that the model distinguished between control and BPD samples (AUC = 0.864, Figure 3F). Additionally, ROC curves were used to assess the association between the individual genes and BPD. The AUC was higher than 0.6 for each gene (Figure 3G). These findings suggest the logistic regression model had a higher level of specificity and accuracy than the individual marker genes.

A diagnostic nomogram based on the five marker PRGs was developed (Figure 4A). The nomogram’s predictive accuracy was highlighted using calibration curves (Figure 4B). The line in the DCA curve remained above the gray and black lines from 0 to 1, suggesting that decisions based on the nomogram model may favor patients with BPD (Figure 4C). The clinical impact curve indicated that the nomogram had a remarkable predictive ability (Figure 4D).

Single-gene GSEA-KEGG pathway analysis was carried out to assess the correlation between BPD and the marker genes. Figures 5A–E presents the six leading pathways that were enriched for each of the marker genes. CHMP7 was closely related to complement and coagulation cascades, oxidative phosphorylation, primary immunodeficiency, ribosome, and splicesome. NLRC4 was related to the cell intestinal immune network for IgA production, complement and coagulation cascades, antigen processing and presentation, primary immunodeficiency, adhesion molecules cams, and ribosome. NLRP2 correlated with the intestinal immune network for IgA production, cell cycle, neuroactive ligand–receptor interaction, olfactory transduction, spliceosome, and valine, leucine, and isoleucine degradation. NLRP6 was related to cytokine–cytokine receptor interaction, proteasome, purine metabolism, ribosome, spliceosome, and valine, leucine, and isoleucine degradation. NLRP9 was related to drug metabolism by cytochrome, maturity-onset diabetes of the young, metabolism of xenobiotics by cytochrome, ribosome, spliceosome, and valine, leucine, and isoleucine degradation.

We then examined the differentially activated pathways between the low- and high-expression groups. Notably, high CHMP7 expression activated glycosphingolipid biosynthesis lacto and neolacto series, complement and coagulation cascades, and folate biosynthesis, whereas the low expression of CHMP7 activated primary immunodeficiency, RNA polymerase, and ribosome. The upregulation of NLRC4 activated the intestinal immune network for IgA production, primary immunodeficiency, and aminoacyl tRNA biosynthesis, whereas the downregulation of NLRC4 activated folate biosynthesis, complement and coagulation cascades, and glycosphingolipid biosynthesis lacto and neolacto series. High NLRP6 expression was associated with protein export, glyoxylate and dicarboxylate metabolism, and butanoate metabolism. Low expression of NLRP6 activated glycosphingolipid biosynthesis lacto and neolacto series, glycosaminoglycan biosynthesis chondroitin sulfate, and complement and coagulation cascades (Figures 6A–C).

The previous findings suggested an association of the immune response with the marker genes. Furthermore, mounting evidence suggests a close association of BPD with the immune microenvironment. Moreover, the CIBERSORT tool was employed to evaluate the association between the immune microenvironment and BPD. The proportions of neutrophils, macrophages M0, and monocytes were elevated in the BPD samples compared to the control samples, whereas the proportions of T cells CD4 naive, T cells CD4 memory resting, T cells CD8, B cells naïve, macrophages M2, and dendritic cells activated were elevated in the control samples compared to the BPD samples (Figure 7A). Furthermore, Pearson’s correlation analysis indicated remarkable negative and positive correlations between NLRP6 and mast cells resting, monocytes, and neutrophils. Strong negative and positive correlations between NLRC4 and macrophages M0, neutrophils, T cells CD4 naive, T cells CD8, T cells CD4 memory resting, and B cells naive were also observed. Additionally, strong positive and negative correlations between CHMP7 and macrophages M0, neutrophils, T cells CD4 naive, T cells CD8, and T cells CD4 memory resting were revealed (Figure 7B).

A ceRNA network based on the five marker genes was constructed. A total of 113 nodes (3 marker genes, 48 miRNAs, and 62 lncRNAs) and 118 edges were present in the network (Figure 8). Five lncRNAs regulated the expression of CHMP7 by binding competitively to hsa-miR-146a-3p. However, hsa-miR-185-5p bound to 10 lncRNAs to exhibit its regulatory role on CHMP7. CTA-414D7.1 regulated NLRP2 expression by binding competitively to hsa-miR-194-5p. LINC00662 regulated NLRP9 expression competitively binding to hsa-miR-708-3p.