From the GEO dataset, we retrieved RNA-seq data (GSE157103) of 126 ARDS patients, namely 100 COVID-19 patients and 26 non-COVID-19 patients, in addition to their clinical data, including gender, age, underlying disease status (diabetes), coagulation (D-dimer, ferritin, CRP, procalcitonin, fibrinogen), and hospital-free days post 45-day follow-up (HFD45), among other parameters. More specifically, HFD45 is defined as the number of days patients lived outside of a hospital from enrollment through death or the end of follow-up (44). The higher HDF45 is, the milder the disease and the better the prognosis will be.

The label-free quantification (LFQ) intensities of 736 proteins of 126 ARDS patients were collected from Mass Spectrometry Interactive Virtual Environment (MassIVE) (MSV000085703). After calculating the logarithm of the FQL intensities, we used the R package “limma” to calculate the log2(Log2Fold of change) (log2FC) of 736 proteins between two PYRclusters. The 736 proteins were sorted from large to small by log2FC (not absolute value). Then, using the “org.Hs.eg.db” and “clusterProfiler” packages, gene set enrichment analysis (GSEA) was performed based on MSigDB gene sets C2, C5 and C7. Significant gene sets were identified when |Normalized Enrichment Score (NES)|>1 and False Discovery Rate (FDR) < 0.25.

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed using the “org.Hs.eg.db” and “clusterProfiler” packages. GO terms and KEGG terms were identified as significantly enriched when p.adjust < 0.05.

Single-sample gene set enrichment analysis (ssGSEA) and “Cibersort” were used to analyze the immune cell infiltration fractions. The former was based on the list of Pan-cancer Immune Metagenes (45, 46).

Based on RNA-seq data of 35 pyroptosis-related genes including AIM2, CASP1, CASP3, CASP4, CASP5, CASP6, CASP8, CASP9, ELANE, GPX4, GSDMA, GSDMB, GSDMC, GSDMD, GSDME, GZMA, GZMB, IL18, IL1B, IL6, NLRC4, NLRP1, NLRP2, NLRP3, NLRP6, NLRP7, NOD1, NOD2, PJVK, PLCG1, PRKACA, PYCARD, SCAF11, TIRAP, and TNF, we divided 100 COVID-19 patients into two clusters (PYRcluster) using the “nmf” package. We determine the k value based on the consensusMap function.

The differentially expressed genes (DEGs) of two PYRclusters were identified using the “limma” package. Using the median of HFD45 as the cutoff value, COVID-19 patients were divided into two groups. The DEGs in these two groups were calculated using the limma package and subsequently applied for PCA analysis. PC1 and PC2 were used to construct the PYRscore (47).

Based on the log2(TPM) of the 570 DEGs between PYRclusters and the HDF45 of each patient, we used the “glmnet” package to build a PYRsafescore model by LASSO regression. We determine the signatures of the model by selecting the lambda value with the smallest mean-squared error by 20-fold-cross-validation. The coefficients of the final signatures were used to calculate the PYRsafescore as follows: protective score = ∑ Coefficient_i × Expression level of signature_i. Using the “caret” package, 100 patients were randomly divided into a training group and a test group with a ratio of 2:1. The model built with the training group data was validated in the test group. We used “ROCR” packages to plot receiver operating characteristic (ROC) curves and to calculate area under the curve (AUC) scores to evaluate model performance.

First, based on the transcription factor targets (TFT) and their gene sets in MSigDB (48), we used the “clusterProfiler” package to enrich the transcription factors from DEGs between PYRclusters.

We then used the “CoRegNet” package to enrich the transcription factor co-regulatory network in COVID-19 patients from PYRcluster1 and PYRcluster2 based on the dataset of transcription factor targets (CHEA, ENCODE, and JASPAR Predicted Transcription Factor Targets, MotifMap and TRANSFAC Predicted Transcription Factor Targets, and TRANSFAC Curated Transcription Factor Targets) from Harmonizome (49, 50). Lastly, the network graph was plotted using the “ggraph” package.

We used the transcription factors that we screened as keywords to search for the corresponding compounds on ChEMBL (51), thus screening active or repressed compounds based on their role in pyroptosis.

The Wilcoxon sum-rank test and the t-test were used to compare different groups, and the Pearson’s product-moment correlation test was used for correlation analysis. All statistical tests were two-sided, and a significant difference was defined as a p-value of 0.05. Power calculations were performed using the following R packages: “pwr” and “rstatix” at sig.level=0.05.

A GEO dataset provided RNA-seq data and clinical data from 126 samples of 100 patients with COVID-19 and 26 patients without COVID-19 (GSE157103) (44). Initially, we assessed the expression levels of pyroptosis-related gene sets of all patients based on prior studies (52).

The expression of the pyroptosis-related genes AIM2, CASP1, CASP3, CASP6, CASP8, CASP9, GSDMA, GSDMC, GZMB, IL6, NLRP3, NLRP7, NOD1, NOD2, SCAF11, and TIRAP was significantly higher in blood leukocytes of COVID-19 patients, indicating that the level of pyroptosis was significantly increased ( Figure 1A ). NLRP3, NLRP7, NOD1, NOD2 are closely related to caspase activation (53, 54).

Subsequently, we outlined the correlation patterns of 35 pyroptosis-related genes in COVID-19 patients to investigate relationships between different pyroptosis-related genes. Although most of the 35 genes have a substantial positive association with other genes, several genes are nevertheless negatively correlated with other pyroptosis-related genes. ELANE, for example, is negatively correlated with AIM2, CASP1, CASP3, CASP4, CASP5, CASP6, CASP8, CASP9, GSDMB, GSDMC, GZMA, IL1B, IL6, NLRP1, NLRP2, NLRP3, NOD1, NOD2, PJVK, PLCG1, PRKACA, SCAF11, TIRAP, and TNF. Therefore, the expression of these pyroptosis-related genes in COVID-19 patients presents complex patterns ( Figure 1B ). Furthermore, each caspase is strongly and positively correlated with the others. For example, the linkage between GSDMA and caspases suggests that its cleavage is related to caspase-3/-6/-8/-9. CASP3 expression is highly linked to GSDME, in line with previous reports on caspase-3 cleavage of GSDME, releasing its activity (22).

By gene set variation analysis (GSVA), we studied changes in the biological function of leukocytes between the two types of patients. Mismatch repair, homologous recombination, replication, cell cycle, and p53 signaling are more enriched in leukocytes of COVID-19 patients, indicating severe cell damage during viral infection and ongoing damage repair ( Figure 2A and Supplementary Table 1 ).

By single-sample gene set enrichment analysis (ssGSEA), we also compared the proportions of 28 immune cell types between the two groups and found that the numbers of various immune cells are much higher in COVID-19 patients than in non-COVID-19 patients, indicating a highly active immune response in COVID-19 patients (46) ( Figure 2B and Supplementary Table 2 ). Surprisingly, the numbers of some immune cells, such as macrophage and cd56^dim natural killer cells, were lower in COVID-19 patients than in the control group. This difference could be due to cell death caused by massive viral infection and increased pyroptosis (17, 55). Human dendritic cells and T lymphocytes (including CD4⁺ and CD8⁺) undergo pyroptosis via the AIM2-Caspase1-gasdermin D and the CARD8-Caspase1-gasdermin D axes (56, 57), respectively. Pyroptosis has also been identified in macrophages and neutrophils (58, 59). Although NK cells have not been documented to undergo pyroptosis on their own, they can participate in this process, playing a key role (60). These results suggests that blood leukocytes either directly undergo pyroptosis or have a synergistic role with pyroptosis in COVID-19 patients.

In summary, blood leukocytes exhibit substantial pyroptotic characteristics in COVID-19 patients.

By non-negative matrix factorization based on 35 pyroptosis-related genes, we clustered the 100 COVID-19 patients into two clusters, PYRcluster1 and PYRcluster2 ( Figure 3A and Supplementary Figure 1A ). The best clustering result was found when k = 2, with no differences in age, gender, days admitted before enrollment, replacement therapy (pre-enrollment) or underlying disorders between the two clusters, which we called PYRclusters ( Supplementary Figures 1B–E , 3A, B ).

We first explored differences in the expression of pyroptosis-related genes between the two PYRclusters ( Figures 3B,C and Supplementary Table 3 ) and found that PYRcluster1 has a higher expression of NLRC4, NLRP6, CASP5, CASP8, CASP9, GSDMC, and AIM2, whereas PYRcluster2 has a higher expression of GPX4, GSDMD, PYCARD, TNF, IL6, NLRP2, NLRP7, GSDMA, CASP6, GSDMB, NOD1, GZMA, and GZMB.

Subsequently, we found that PYRcluster2 has a higher hospital-free days post 45-day follow-up (HFD45), ventilator-free days and a lower proportion of mechanical ventilation than PYRcluster1 ( Figures 3D, E and Supplementary Figure 3C ), which suggests a better prognosis. The blood D-dimer level of PYRcluster1 is significantly higher than that of PYRcluster2 ( Figure 3F ), indicating hypercoagulation. Although albumin and hemoglobin of the two PYRclusters are mostly below the normal range (green dashed line), PYRcluster1 deviates further from the normal range than PYRcluster2. Other clinical features do not differ significantly between the two clusters ( Supplementary Figures 2A–C ). Based on these results, the highly expressed pyroptosis-related genes of PYRcluster1 may be associated with a poor prognosis. In line with our results, AIM2 and NLRC4 deficiency in dogs and cats provide a protective effect against SARS-CoV-2 by reducing cytokine storm-induced host damage (43).

We further determined differentially expressed genes (DEGs) using the “limma” package, with cutoff criteria of |logFC| >1 and p=0.05, totaling 570 DEGs, and employed Gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichments ( Figure 4A and Supplementary Figure 2E ). In addition, we determined the logFC of 736 proteins between two PYRclusters (FC=PYRcluster2/PYRcluster1) using “limma” from the proteomic data from the same study as the transcriptome data. Subsequently, we employed gene set enrichment analysis (GSEA) based on the log2FC of these proteins ( Supplementary Figure 2D , Supplementary Table 4 ).

Immune responses such as antigen recognition/presentation, immune cell activation, migration, and replication are relatively enhanced in PYRcluster2 ( Figure 4A and Supplementary Figure 2E ). Both transcriptome and proteomic results revealed hyperinflammation in PYRcluster1 ( Supplementary Figure 2D , Supplementary Table 6 ). KEGG pathway enrichment data showed that “Coronavirus disease-COVID-19” pathways are markedly upregulated in PYRcluster2, underscoring a highly activated immune process unique to PYRcluster2. PYRcluster1 has a higher NO synthesis level than PYRcluster2, which may be related to the antiviral capacity of its patients (61–63). The highly activated stress response pathway may associated with severe damage caused by the virus and lead to higher levels of immune cell apoptosis in PYRcluster1 ( Supplementary Figure 2D , Supplementary Table 6 ). Furthermore, PYRcluster1 has a markedly increased expression of cytokines, including IL-1, IL-6, IL-8, IL-10, and TNF ( Supplementary Table 6 ), which are typical components of a cytokine storm and result in a poor prognosis (64, 65). Consistent with an increased D-dimer ( Figure 3F ), the neutrophil extracellular traps (NETs) that promote blood coagulation are highly expressed in PYRcluster1, promoting venous thrombosis and leading to poor prognosis (66, 67). In addition, AIM2, CXCR2 and UBE2W, which were significantly highly expressed in PYRcluster1, were included in the 20 genes associated with distinctly methylated CpG sites between mild and severe COVID-19 patients (68). Their odds ratios were all greater than 1, indicating that their downregulation is beneficial to COVID-19 patients. In conclusion, the two PYRclusters of COVID-19 patients exhibited distinct pyroptotic patterns and clinical features.

SARS-CoV-2 viruses infect leukocytes and lead to immunodeficiency (69–71). We speculate that leukocyte pyroptosis helps to destroy the virus protective niche and release viruses from cells, thereby enhancing viral clearance and immune recovery. As a result, viruses released by pyroptosis may be further removed by phagocytic cells via phagocytosis.

We first assessed the proportion of different immune cells in the two PYRclusters by ssGSEA. Immune cells highly associated with antivirals, such as activated B, CD4/8⁺ T, Treg, and NK cells are highly expressed in PYRcluster2 ( Figure 4B and Supplementary Figure 4A , Supplementary Table 7 ). PYRcluster1 had a higher proportion of pro-inflammatory neutrophils. Using “Cibersort”, we discovered PYRcluster2 had a higher proportion of anti-inflammatory M2 macrophages ( Supplementary Figure 4A ). The hemogram percentages of several cells, such as neutrophils, lymphocytes and monocytes, were consistent with the results of ssGSEA ( Supplementary Figure 3F ). We then used the “estimate” package to score the two clusters and found that PYRcluster2 shows a greater increase in immune cell infiltration than PYRcluster1 ( Figure 4C ). As shown in Figure 4D , the expression of characteristic pyroptosis-related genes of PYRcluster1 were significantly positively correlated with its expression of characteristic immune cells. For example, PYCARD, GPX4, GSDMD, GZMA, TNF, NOD1, IL6, NLRP7, CASP6, GSDME, PJVK, GSDMA, and NLPR2 are positively correlated with NK, NKT, regulatory T, γδ T, activated B, and Th17 cells, MDSCs, and monocytes ( Figure 4D ). Both these pyroptosis-related genes and immune cells are highly expressed in PYRcluster2 ( Figures 3C , 4B ). Furthermore, leukocytes have stronger phagocytic activity in GO enrichment results in PYRcluster2 than in PYRcluster1 ( Supplementary Table 6 ). In contrast, pyroptosis in PYRcluster1 may produce several pathogen- and damage-related molecular patterns that increase cytokine storm, leading to multiple organ failure and poor prognosis.

We subsequently performed unsupervised clustering of all COVID-19 patients using 570 DEGs and all genes, respectively, yielding two more clusters: DEG and All-Gene clusters. The heatmap demonstrates that these additional clusters match a previous clustering based on the 35 pyroptosis-related genes ( Figure 5A and Supplementary Figure 4B–D ). These data suggest that distinct patterns of pyroptosis occur in COVID-19 patients, which can be represented by 35 pyroptosis-related genes.

To better elucidate differences in pyroptotic patterns between PYRcluster1 and PYRcluster2 and their correlation with prognosis, we created a pyroptosis score (pyrscore) ( Supplementary Figure 4E ). As shown in Figure 5B , PYRcluster1 has a higher score than PYRcluster2. HFD45 and ventilator-free days are negatively correlated with pyroptosis scores ( Figures 5C, D ), whereas sofa, APACHE-II, D-dimer, and CRP levels are positively correlated ( Supplementary Figures 5A–D ). In conclusion, PYRcluster1 and PYRcluster2 have different levels of immune response to SARS-COV-2 and pyroptotic patterns, resulting in distinct prognoses.

Given that different pyroptotic patterns may have a significant impact on the prognosis of COVID-19 patients, we created a PYRsafescore model based on the HFD45 of COVID-19 patients and DEGs across PYRclusters by least absolute shrinkage and selection operator (LASSO) regression analysis. In a ratio of 2:1, 100 COVID-19 patients were divided into a training group and a test group, and the model was obtained in the training group. “Lambda-min” was chosen as the best value in the cross-validation procedure ( Figure 5E and Supplementary Figures 5E ). Lastly, based on the log2 value of the expression level of 10 genes, we established the following scoring model:

The area under the ROC curve (AUC) values of the model in the training and test groups were 0.907 and 0.879, respectively, indicating that our predictive model performs well ( Figure 5F ). HFD45 and PYRsafescore are positively correlated ( Figure 6A ), which suggests that a higher PYRsafescore indicates a better prognosis. In addition to HDF45, the correlation of other clinical variables with PYRsafescore, including sofa, ventilator-free days, APACHE-II, and CRP levels, also demonstrates that our model works effectively ( Figures 6B, C and Supplementary Figures 6A, B ). Additionally, the expression of these 10 genes between the two PYRclusters also has significant differences ( Figure 6D ). KEGG analysis indicates that they were closely related to the patient’s immune response ( Figure 6E ).

In conclusion, our results proved that the newly created prognosis model has a considerable clinical predictive value.

Using the “clusterProfiler” package, we first enriched DEGs across PYRclusters for transcription factors based on MSigDB Collections: “regulatory target gene sets” ( Figure 6F ). We further enriched the transcription factor regulatory networks for all transcriptome data of the two PYRclusters using the “CoRegNet” package. PYRcluster1 has the regulatory network with MNDA, TSC22D3, HMGB2, FOS, EEF1A1, TRIM22, NFKBIA as transcription factors, and PYRcluster2 has the regulatory network with FOS, MNDA, EEF1A1, DAZAP2, DDIT3, HCLS1, NFKBIA, TSC22D3, PTMA, and TRIM22 ( Figure 7A and Supplementary Figure 6C , Supplementary Table 8 ). Notably, both PYRcluster1 and PYRcluster2 are regulated by FOS, EEF1A1, MNDA, and TRIM22.

Once SARS-COV-2 infects cells, EEF1A1 is critical for viral replication. Drug targeting EEF1A has robust antiviral effects in vitro (72). The nuclear factor-kappaB (NF-kappaB)/REL family of transcription factors are activated in response to DNA damage to regulate inflammation and apoptosis resistance (73, 74). Since NFKBIA is highly expressed in COVID-19 patients (75) and tripartite motif-containing (TRIM) 22 can activate NF-kappaB (76) to protect the host from viral infection (77), these two genes have a huge impact on immune disparities between the two PYRclusters. MNDA is a member of the family of hematopoietic interferon (IFN)-inducible nuclear proteins that promotes the degradation of the anti-apoptotic factor MCL-1 and apoptosis in myeloid cells (78, 79).

We then observed the correlation between these transcription factors and the expression of differentially expressed pyroptosis-related genes between PYRclusters and the clinical characteristics of all patients. The correlations between transcription factors and pyroptosis-related genes vary greatly, depending on the PYRclusters ( Figure 7B ). Pyroptosis-related genes highly expressed in PYRcluster1 are positively correlated with FOS, MAML1, DAZAP2, STAT6, STAT5B, ETS2, HCLS1, TSC22D3, HMGB2, MNDA, TRIM22, NFKBIA, and DDIT3, whereas genes highly expressed in PYRcluster2 are positively correlated with ZNF623, PSMB5, PTMA, EEF1A1, ETS1, CEBPA, and USF2. Moreover, some transcription factors positively correlated with pyroptosis-related genes in PYRcluster2 also have a positive effect on prognosis ( Figure 7C ).

Taken together, these findings show that distinct pyroptotic patterns may result from different upstream transcriptional regulation pathways. In light of this, we searched the ChEMBL database for compounds that promote or inhibit appropriate transcription factors based on diverse roles in prognosis (51). We screened a total of 112 compounds ( Supplementary Table 9 ) and found that CHEMBL348436 (also known as Cirsimaritin) has the potential to regulate blood leukocyte pyroptosis in COVID-19 patients while simultaneously improving prognosis through FOS and NFKB1A inhibition. In fact, drugs targeting FOS have therapeutic effects in COVID-19 patients (80).