Two peripheral-blood gene expression datasets (GSE65682 and GSE95233) and their corresponding clinical data were obtained from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/gds/). The GSE65682 dataset and GSE95233 dataset were based on GPL13667 and GPL570 platforms, respectively. The GSE65682 dataset comprised 479 sepsis patients with complete data for survival status within 28 days and 42 healthy controls. Sepsis patients in GSE65682 were separated into the discovery cohort (n = 263) and validation cohort (n = 216) for its original investigation and we used these two cohorts for survival analysis in our study. The GSE95233 dataset comprised 51 sepsis patients and 22 healthy controls. Only the gene expression data of the blood sample collected at admission to the intensive care unit was used for analysis in the present study.

The gene probe was transformed into gene symbol by using the corresponding annotation profile in each dataset. We used ‘limma’ package in R software to quartile normalized for all gene expression values and generate normally distributed expression values. For multiple same probes, the finial gene expression value was determined by calculating the average expression value.

A total of 60 pyroptosis-related genes were identified from gene set enrichment analysis (GSEA) website (http://www.gsea-msigdb.org/gsea/index.jsp) and previous literature (17–21) ( Supplementary Table 1 ). Differentially expressed genes (DEGs) between sepsis and healthy samples in GSE65682 were screened using the ‘limma’ package, with log₂|fold change| > 0.5 and adjust P value< 0.05 set as cut-off criteria to screen DEGs (22, 23). Pyroptosis-related DEGs were identified by intersecting the DEGs with the pyroptosis-related genes. A protein-protein interaction (PPI) network analysis was conducted by using the STRING website (https://cn.string-db.org/) to further explore the interaction between these pyroptosis-related genes. Based on the expression value of pyroptosis-related DEGs, we used “ConsensusClusterPlus” package to identify the molecular subtype of sepsis. The pam algorithm with euclidean distance was used, and the samples were iterated 1000 times. The k value was increased from 2 to 6 to identify the optimal clusters.

Univariate Cox regression analysis was performed in the discovery cohort to evaluate the prognostic value of each pyroptosis-related DEGs. To avoid omissions, we set a P-value lower than 0.2 as a significant cut-off value (24). Pyroptosis-related DEGs with a significant correlation to survival status were selected as candidate genes for further investigation. Then, we used LASSO-Cox regression via ‘glmnet’ packages to screen the candidate genes and construct the prediction model. The penalty coefficient λ was determined by using the minimal criteria and candidate genes with a regression coefficient unequal to zero were included in the final model. The risk score of each case was calculated according to the following formula: Risk score  = ∑ i k Coefficient i  ×  Gene Expression value of i   . After that, sepsis patients were divided into low- and high-risk groups based on the median risk score, and the survival time between the two risk groups was compared using Kaplan-Meier analysis. Time-dependent receiver operating characteristic (ROC) curve analysis was also performed by using ‘survival’, ‘surviminer’, ‘timeROC’ packages to assess the discrimination ability of the risk score. The risk score was further validated in the validation cohort and the whole sepsis patients in GSE65682 (test cohort). Furthermore, we compared the 28-day mortality rate of sepsis patients from the GSE95233 dataset who were classified as low- or high-risk based on risk score (GSE95233 cohort). We were unable to perform the Kaplan-Meier analysis because the GSE95233 dataset lacked “time to event” data. A ROC curve was depicted via ‘pROC’ package.

Clinical information (age and sex) was extracted from patients in the GSE65682 dataset. These variables and risk score were analyzed together in univariate and multivariate Cox regression analyses. The ‘rms’ package was used to create a nomogram based on independent variables for visualization and potential clinical use in predicting the prognosis of sepsis patients. The ROC curve and calibration curve were used to assess the nomogram’s performance.

We used ‘limma’ package again and based on the same criteria (log₂|fold change| > 0.5 and adjust P value< 0.05) to screen DEGs between low- and high-risk groups in the sepsis patients from GSE65682. To explore the DEGs-related signal pathways and biological function, “clusterprofiler” package (25) was applied to perform the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analysis. To quantify the relative proportion of immune cell infiltration, we used the CIBERSORT algorithm (26) with 1000 permutations to calculate 22 types of immune cell composition for each sample. The composition of these immune cells between low- and high-risk groups was compared via wilcoxson test. In addition, we performed a correlation analysis between the differentially immune cells and the pyroptosis-related genes.

We further investigated whether the selected prognostic pyroptosis-related genes also have potential value in the diagnosis of sepsis. The performance of these genes was evaluated using the ROC curves in the GSE65682 dataset and GSE95233 dataset.

All statistical analyses were performed using R software (version 4.1.0) and R studio (Version 1.2.5042). The wilcoxon test was used to compare the gene expression level between sepsis patients and healthy individuals, and the composition of immune cells between low- and high-risk groups. Pearson chi-square test was applied to compare categorical variables. LASSO-Cox regression was used for candidate genes selection. The Kaplan-Meier method and log-rank test were used to compare the survival rate between the low- and high-risk groups. Univariate and multivariate cox regression analyses were used to assess the independent prognostic variables. A two-tailed P value<0.05 was considered statistically significant except for a certain P value was set.

A total of 3469 DEGs were identified from the GSE65682 dataset between sepsis and healthy samples, including 1571 upregulated genes and 1898 downregulated genes ( Figure 1A ; Supplementary Table 2 ). After intersecting with the pyroptosis-related genes, 16 pyroptosis-related DEGs were obtained ( Figure 1B ). Among them, the expression level of 7 genes (MYD88, NLRP3, TLR2, CASP5, NLRC4, ELANE, AIM2) were upregulated, while the expression level of 9 genes (GZMB, CHMP7, NLRP1, IRF1, PLCG1, SCAF11, AKT1, GSDMB, IRF2) were downregulated ( Figure 1C ). The result of the PPI analysis was presented in Figure 1D . There were 40 interaction relationships between these pyroptosis-related DEGs.

According to the empirical CDF value, k = 2 was found to be the most acceptable point for the consensus cluster with the most distinct differences between clusters ( Figure 2A ; Supplementary Figure 1 ). Some of the expression levels of the pyroptosis-related DEGs were different between the two clusters ( Figures 2B, D ), and sepsis patients in cluster 1 had a worse prognosis than patients in cluster 2 (P = 0.0099) ( Figure 2C ).

Based on the univariate cox regression analysis, 10 out of 16 genes (GZMB, CHMP7, NLRP1, IRF1, SCAF11, IRF2, MYD88, CASP5, ELANE, AIM2) met P< 0.2 and were selected as prognostic candidate genes ( Figure 3A ). Of the 10 prognostic candidate genes, 1 gene (ELANE) was associated with increased risk with (HR > 1), while the remaining 9 genes (GZMB, CHMP7, NLRP1, IRF1, SCAF11, IRF2, MYD88, CASP5) were associated with lower risk (HR< 1). By performing LASSO-Cox regression analysis with the above candidate genes, a subset of six genes (GZMB, CHMP7, NLPR1, MYD88, ELANE, AIM2) were determined to develop a prognostic risk score based on the minimal criteria of λ ( Figures 3B, C ). The calculation of the risk score for each sample was according to the formula as follows: Risk score = [(-0.01470407 x GZMB expression value) + (-0.63970285 x CHMP7 expression value) + (-0.12447921 x NLRP1 expression value) + (-0.18172740 x MYD88 expression value) + (0.07255413 x ELANE expression value) + (-0.16451424 x AIM2 expression value)]. After calculating the median risk score, 263 patients were stratified into two risk groups (131 in the low-risk group and 132 in the high-risk group) ( Figure 3D ). Patients in the high-risk group had more death ( Figure 3E ) and the Kaplan-Meier curve showed that patients in the low-risk group had a higher survival rate than those in the high-risk group (P = 0.0011, Figure 3F ). Time-dependent ROC analysis of the risk score revealed that the area under the curve (AUC) was 0.70 (95% CI: 0.60 to 0.80), 0.68 (95% CI: 0.59 to 0.76), and 0.65 (95% CI: 0.57 to 0.73) for 7-, 14-, and 28-day survival, respectively ( Figure 3G ). Both univariate and multivariate cox regression analysis indicated that the risk score was an independent predictor for prognosis in sepsis patients (Univariate: HR = 3.78, 95% CI: 2.33 to 6.4, P< 0.001; Multivariate: HR = 3.72, 95% CI: 2.16 to 6.39, P< 0.001; Figure 3H ).

Patients in the validation cohort and the test cohort were divided into two groups based on the median risk score, respectively ( Figures 4A, E ). Patients in the high-risk groups also occurred more death events ( Figures 4B, F ) and Kaplan-Meier curve showed that patients in the low-risk group had significantly higher survival rates than those in the high-risk group (In the validation cohort: P = 0.013, Figure 4C ; In the test cohort: P< 0.0001, Figure 4G ). AUC of the time-dependent ROC curve was 0.66 (95% CI: 0.53 to 0.79), 0.63 (95% CI: 0.52 to 0.73), and 0.64 (95% CI: 0.55 to 0.72) for 7-, 14-, and 28-day survival in the validation cohort ( Figure 4D ) and 0.69 (95% CI: 0.61 to 0.77), 0.65 (95% CI: 0.59 to 0.72), and 0.64 (95% CI: 0.59 to 0.70) for 7-, 14-, and 28-day survival in the test cohort ( Figure 4H ). According to the result of univariate and multivariate cox regression, the risk score also could be an independent prognostic factor in sepsis patients in these two cohorts ( Supplementary Figure 2 ). In addition, patients in the low-risk group also had a higher survival rate than those in the high-risk group in the GSE95233 cohort (P = 0.047, Figure 4I ). The AUC under the ROC curve was 0.687 (95% CI: 0.531 to 0.842) ( Figure 4J ).

In order to potentially clinically used the risk score and predicted more precisely the prognosis of sepsis patients. We created a nomogram using the patient’s age (P = 0.029, Supplementary Figure 2B ) and risk score. ( Figure 5A ). The AUC of the nomogram for predicting 7-, 14-, and 28-day of survival were 0.69 (95% CI: 0.61 to 0.77), 0.69 (95% CI: 0.62 to 0.75), and 0.67 (95% CI: 0.61 to 0.72) ( Figure 5B ). The calibration curve indicated a good calibration between predicting probability and actual probability ( Figures 5C–E ).

A total of 485 DEGs were identified between the low- and high-risk groups. Among them, 164 genes were downregulated and 321 genes were upregulated ( Supplementary Table 3 ). On the basis of these DEGs, KEGG and GO analyses were performed. The results showed that the DEGs were mainly correlated with NOD-like receptor signaling pathway, hematopoietic cell lineage, and response to the virus ( Figures 6A, B ).

Based on the results of the functional analysis. We further compared immune infiltration between low- and high-risk groups by using the CIBERSORT algorithm. In the low-risk group, neutrophils, B cells memory, and mast cells activated were significantly higher than that in the high-risk group, while T cells CD4 naive, macrophages M0, mast cells resting, eosinophils, macrophages M2, plasma cells, and T cells CD4 memory resting were significantly lower ( Figure 7A ). The correlation analysis revealed that pyroptosis-related genes were significantly associated with many immune cells. All six genes were correlated with neutrophils, with AIM2 (r = 0.34), MYD88 (r = 0.4), and NLRP1 (r = 0.47) showing a positive correlation, and CHMP7 (r = -0.2), ELANE (r = -0.22), and GZMB (r = -0.17) showing a negative correlation ( Figure 7B ).

According to the ROC curves, four out of six genes (GZMB, CHMP7, NLRP1, and AIM2) had good diagnostic value in the diagnosis of sepsis, with AUC > 0.9 in both the GSE65682 and GSE95233 datasets ( Figures 8A, B ). In both datasets, the expression level of AIM2, ELANE, and MYD88 were significantly higher in sepsis patients compared to healthy individuals, while the expression level of CHMP7, GZMB, and NLRP1 were significantly lower. ( Figures 8C, D ).