Machine Learning Screening and Validation of PANoptosis-Related Gene Signatures in Sepsis

Introduction

Sepsis, a syndrome characterized by life-threatening organ dysfunction and a dysregulated host immune response to infection,^1,2 poses a significant global health concern. According to a recent global burden of disease report, there were 48.9 million sepsis cases, with a mortality rate of 22.5%, accounting for nearly 20% of all global deaths.³ The substantial prevalence and mortality associated with sepsis, coupled with the absence of effective diagnostic and treatment approaches, present a challenge for proper treatment to intensive care physicians.⁴ While a growing number of biomarkers have been associated with sepsis, only few have been thoroughly evaluated.^5,6 Hence, the accurate molecular-level identification of sepsis and the development of multivariate predictive models hold significant clinical significance.

Programmed cell death (PCD) represents a genetically determined mechanism of active cell death for eradicating pathogens and preserving cellular equilibrium. It encompasses various forms, such as pyroptosis, apoptosis, and necroptosis.^7,8 Recognizing that targeting individual cell death pathways in isolation may not yield the desired therapeutic outcomes, recent research has increasingly focused on understanding the redundancy and interplay among these pathways. A recent conceptual development is PANoptosis, which emphasises the interconnectedness and cooperation among pyroptosis, apoptosis, and necroptosis, triggered by specific stimuli and regulated through the PANoptosome complex.⁹ Recent investigations have linked PANoptosis activation to responses triggered by viruses, bacteria, and fungi, as well as its association with autoimmune diseases, cytokine storms, and cancer.¹⁰ Furthermore, several studies have indicated the pivotal role of PANoptosis in multiorgan dysfunction resulting from sepsis. Inhibition of PANoptosis has shown promise in mitigating acute lung injury¹¹ and sepsis-associated encephalopathy.¹² However, the specific role of PANoptosis in sepsis is still unknown and warrants further exploration.

In our study, we carried out the first systematic exploration of differentially expressed PANoptosis-related genes (PRGs) and immune infiltration in sepsis. Depending on the landscape of 19 PRGs expression, we categorized sepsis samples into two PANoptosis-associated clusters and analyzed their immune cell differences. We then identified PANoptosis-specific DEGs and established a predictive model using various machine learning algorithms to distinguish patients based on their molecular clusters. The model’s performance was validated through an external dataset, nomograms, survival analysis, quantitative real-time PCR, and western blot. Collectively, this study identified promising diagnostic biomarkers through machine learning, offering a novel perspective for potential promising developments to treat sepsis more effectively.

Materials and Methods

Data Collection and Processing

We obtained two sepsis gene expression datasets from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/gds/): GSE65682 (from the GPL13667 platform) and GSE95233 (from the GPL570 platform). GSE65682 was used as the training dataset, and we focused on 479 sepsis samples with complete survival data and 42 healthy controls in order to analyze (Supplementary Tables 1 and 2; Supplementary Material 2).¹³ GSE95233 was used for validation, and the gene expression data from blood samples collected at day 1–2 after sepsis in the ICU were extracted for analysis.¹⁴ Nineteen PRGs were retrieved from prior studies.^15–19

Immune Cell Infiltration

We used CIBERSORT (https://cibersortx.stanford.edu/), a gene-based deconvolution algorithm, to quantify relative scores for 22 immune cells using standardized gene expression data.²⁰ The correlation coefficients between PRGs and immune cell properties associated with sepsis were visualized using the “corrplot” R package.

Consensus Clustering Analysis

We utilized the “ConsensusClusterPlus” R package to perform clustering analysis based on PRGs expression and classify sepsis samples into two subgroups, while the “stats” R package was used for PCA. One hundred and nine cluster-related DEGs in all were identified comparing the cluster2 samples to the cluster1 samples by the “limma” R package (p-value < 0.05, | log₂Fold Change | > 1).

Functional Enrichment Analysis

Using the “clusterprofiler” R package, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were conducted to clarify the potential biological functions and signaling pathways involved in these DEGs.

Weight Gene Co-Expression Analysis (WGCNA)

Weight gene co-expression analysis was performed to identify modules associated with PRGs in sepsis by using the “WGCNA” R package.²¹ The approach involved several steps: first, we performed sample clustering to eliminate outliers using the “hclust” function. Second, the “pickSoftThreshold” function was used to establish the required soft thresholding power (β) for achieving a scale-free topology, and we transformed the weighted adjacency matrix into a Topological Overlap Matrix (TOM). Next, the “pickSoftThreshold” and “plotDendroAndColors” functions were used to identify and visualize modules. The clustering dendrogram of genes, displaying the dissimilarity of topological overlaps, was generated, with each module assigned a distinct random color. Module-trait correlation analysis was performed to reveal that the ME red module as the most significant one related to PRGs in sepsis. To investigate further, we identified driver genes in the ME red module based on Gene Significance (GS) and Module Membership (MM) calculations and selected genes with GS > 0.2 and MM > 0.8 for further analysis.

Identification of Differentially Expressed Genes

The “limma” R package was used to screen DEGs between control and sepsis samples in the GSE65682 dataset using the screening criteria |log₂Fold Change | >1 and p-value < 0.05.

Construction of Predictive Model

We constructed four machine learning models, namely the Random Forest model (RF), Support Vector Machine model (SVM), Generalized Linear Model (GLM), and eXtreme Gradient Boosting (XGB) by the “caret” R package, which were then subjected to 10-fold cross-validation, utilizing default parameters. Then the “DALEX” R package was utilized to analyze the residual distribution and feature importance. In addition, we used the “pROC” R package to create visualizations of the Area Under Receiver Operating Characteristic (ROC) curves and assess the classification performance.

Construction of a Nomogram Model and Survival Analysis

To assess five key PANoptosis-related sepsis predictor genes, we construct a nomogram model using the “rms” R package. The calibration curve and decision curve analysis (DCA) were used to calculate the nomogram model’s predictive capability. We performed a Kaplan–Meier (K-M) analysis using the “survival” R package and the p-value of survival curves was calculated by the “survminer” R package.

Cell Culture and Treatment

Human Umbilical Vein Endothelial Cells (HUVECs) were used to establish an in vitro model.²² HUVECs and Dulbecco’s Modified Eagle’s Medium (DMEM) were obtained from Procell Co. Ltd (Wuhan, China). The HUVECs were cultured at 37°C with 5% CO2 using standard methods and treated with 1 μg/mL LPS (Sigma, St. Louis, MO, USA) for 24 hours.

RT-PCR Validation

Total RNA was isolated from cells using Universal RNA Purification Kit (EZBioscience, USA) and reverse transcribed into cDNA by Hifair III 1st Strand cDNA Synthesis Kit (Yeasen, China). Real-time PCR was conducted with Hieff UNICON Universal Blue qPCR SYBR Green Master Mix (Yeasen, China). The primer sequences utilized for PCR amplification were synthesized by Servicebio (Wuhan, China) and presented in Table 1.

Table 1 Sequences of the Primers Used in This Study

Western Blot and ELISA

Total cellular proteins were lysed with RIPA lysis buffer (Beyotime, China). Total proteins were electrophoresed on a 10% SDS-PAGE gel and transferred to a PVDF membrane. The membrane was then closed with 5% skimmed milk and incubated with primary antibodies against DDX60 (dilution 1:1000, Proteintech, China), IFIT3 (dilution 1:2000, Proteintech, China), OAS2 (dilution 1:1000, Proteintech, China), UBE2L6 (dilution 1:1000, Proteintech, China) and GAPDH (dilution 1:50000, Proteintech, China) at 4°C overnight. After washing with TBST, the membrane was incubated with the secondary antibody for one hour at room temperature. Bands were displayed using ChemiDoc XRS (Bio Rad, USA) and optical density was analysed using Image J software. The levels of IL-6 (4A Biotech, China) and TNFa (BOSTER, China) in cell supernatants were measured using enzyme-linked immunosorbent assay (ELISA) kits.

Statistical Analysis

Statistical analysis was performed using R (version 4.3.0) and GraphPad Prism (version 9.0.0) software. One-way ANOVA and t-test were used to compare differences between multiple groups (≥2 groups) and between the two groups, respectively. Statistical significance was set at p<0.05.

Results

Differentially Expressed PRGs in Sepsis and Correlation Analysis Between Genes

Flowchart of our study was shown in Figure 1. The chromosomal locations of copy number variation (CNV) alterations in PRGs are illustrated in Figure 2A. Among them, nine genes, specifically AIM2, NLRP3, PIPK1, GSDMD, FADD, RIPK3, MLKL, PSTPIP2, and TAB3, were upregulated while six genes, namely ZBP1, CASP8, CASP6, PARP1, IRF1, and TRADD were downregulated in sepsis (Figure 2B and C). The Circos plot further shows the interconnectivity between these PRGs with differential expression (Figure 2D).

Figure 1 Flowchart of our study.

Figure 2 Identification of PRGs in sepsis. (A) Locations of CNV variants in PRGs on chromosomes. (B) Boxplot and (C) heatmap demonstrating differences in expression levels of PRGs between healthy controls and sepsis samples. (D) Circos plot displaying the interconnectivity among PRGs. *p < 0.05; ***p ≤ 0.001.

Correlation of PRGs with Immune Infiltration

We conducted immune infiltration analysis to further elucidate potential differences in the immune landscape. Results revealed higher levels of infiltration in sepsis samples, especially plasma cells, T cells, NK cells, and macrophages (Figure 3A), suggesting that alterations in the immune system might play an essential part in sepsis. Additionally, correlation analysis results indicated significant associations between almost all types of infiltrating immune cells and PRGs (Figure 3B). These findings indicate that PRGs may serve as crucial factors influencing the immune infiltration status in sepsis samples.

Figure 3 Correlation of PRGs with immune infiltration in sepsis. (A) Boxplot illustrating the differentiation in immune infiltration. (B) Heatmap displaying significant correlations between PRGs and immune infiltration cells. *p < 0.05; **p ≤ 0.01; ***p ≤ 0.001.

Identification of Two Sepsis Subtypes Based on PANoptosis-Related Genes

Using a consensus clustering algorithm, we successfully identified two distinct subtypes of sepsis to determine PANoptosis-related expression patterns in sepsis. The optimal clustering results were achieved when setting the k value to 2 (Figure 4A and B). PCA showed that there were significantly different among these clusters (Figure 4C). Furthermore, the expression levels of most PRGs were increased in Cluster1 compared to Cluster2, such as AIM2, CASP8, IRF1, RIPK1, GSDMD, RIPK3, MLKL, PSTPIP2, and ZBP1 (Figure 4D and E).

Figure 4 PRGs-based two different clusters of sepsis. (A) Unsupervised consensus clustering matrix and optimal clusters. (B) CDF curves with consensus scores between k = 2–9. (C) PCA based on clustering results. (D) Heatmap and (E) boxplot showing the relationship between two clusters and expression of PRGs in sepsis. **p ≤ 0.01; ***p ≤ 0.001.

Biological Characteristics Between Two Sepsis Clusters

One hundred and nine cluster-related DEGs were defined to analyze PANoptosis-related biological functions, signaling pathways and differences in immune cell infiltration between clusters in sepsis. GO analysis demonstrated that cluster-related DEGs were significantly connected to virus immune- related pathways that were upregulated, such as defense response to virus, negative regulation of viral process, and defense response to symbiont (Figure 5A). KEGG analysis indicated that these genes were involved in immune-related pathways associated with NOD-like receptor signaling, Coronavirus disease (COVID-19), and RIG-I-like receptor signaling pathways (Figure 5B). Immune infiltration analysis showed that Cluster 1 had relatively greater abundances of immunocytes, especially memory B cells, helper follicular T cells, and macrophages (Figure 5C).

Figure 5 Biological characteristics between two clusters. (A) GO enrichment analysis of cluster-related DEGs. (B) KEGG pathway enrichment analysis of cluster-related DEGs. (C) Boxplot illustrating how immune infiltration varies between clusters. *p < 0.05; **p ≤ 0.01; ***p ≤ 0.001.

WGCNA Co-Expression Analysis of Two Sepsis Clusters

Using the WGCNA algorithm, we identified important modules tightly connected to PANoptosis clusters and cluster-related hub module genes in sepsis. Using 0.9 as the scale-free R2 and 9 as the soft power parameter, the samples remained after eliminating the outlier samples were used to identify co-expressed modules (Figure 6A and B). Using the dynamic cutting algorithm, we obtained 9 distinct co-expression modules comprising 2715 genes (Figure 6C and D). These genes were subsequently used to analyze the co-expression similarity and adjacency with cluster-specific clinical characteristics. Lastly, we found that the red module displayed the strongest association with Cluster 2, encompassing 169 genes and that the red module genes exhibited an excellent connection with the selected module (Figure 6E and F).

Figure 6 WGCNA co-expression analysis of two sepsis clusters. (A) Detection of outliers. (B) Determination of the soft power parameter. (C) Tree dendrogram of modules. Different co-expression modules are represented by different colors. (D) Correlation analysis be 9 modules. (E) Heatmap displaying the relationship between phenotypes and module genes. The red module shows significant correlation with PRGs with sepsis (R =0.77, p < 0.001). (F) Scatterplot between module membership (MM) in red module and the gene significance (GS) for Cluster 2.

Construction and Assessment of Machine Learning Models

One thousand two hundred and forty-two DEGs between healthy control samples and sepsis samples are identified in the GSE65682 (Figure 7A and B). As displayed in the Venn diagram of Figure 8A, we identified 12 characteristic PANoptosis-specific DEGs in sepsis, intersecting 169 module hub genes of PANoptosis clusters with 1242 DEGs. Four machine learning models were established, based on these 12 PANoptosis-specific DEGs, and then we assessed the distribution of residuals with each model (Figure 8B and C). Based on the results of the Root Mean Square Error (RMSE) ranking, we obtain the top 10 important characteristic variables of these models (Figure 8D). Subsequently, the ROC curves were used to validate the discriminative performance of each model (Figure 8E). Given its lower residual and higher AUC, the SVM model was considered the best-performing classifier associated with PANoptosis in sepsis. Next, we validated five key PANoptosis-related sepsis predictor genes (DDX60, TAP2, UBE2L6, IFIT3, and OAS2) on an external dataset (GSE95233) and the ROC curve (AUC = 0.989) demonstrates that the model does have satisfactory performance (Figure 8F). These results collectively indicate the efficacy of our diagnostic model in distinguishing between sepsis and normal samples.

Figure 7 Identification of DEGs between control and sepsis samples. (A) Volcano plot of the DEGs. (B) Heatmap showing the differential expression of the top 30 upregulated genes and top 30 down-regulated genes.

Figure 8 Construction and assessment of machine learning models. (A) The intersections between cluster-related hub module genes and DEGs in control and sepsis samples. (B) Reverse cumulative distribution map of RF, SVM, GLM and XGB model residuals. (C) Boxplots showing the residuals of each model. (D) The important features of four machine models. (E) ROC analysis of machine learning models in the testing cohort. (F) ROC analysis of five key PANoptosis-related sepsis predictor genes in the GSE95233.

Construction of the LPS-Induced HUVECs and Verification of Gene Expression

In cell supernatants, we observed an increase in the expression of inflammatory cytokines (IL-6 and TNF-α) with increasing concentrations of LPS stimulation (Figure 9A and B). Therefore, subsequent experiments were performed under 1ug/mL LPS stimulation for 24 hours. The mRNA levels of DDX60, UBE2L6, IFIT3, and OAS2 were all significantly downregulated in LPS-induced HUVECs (Figure 9C). At the protein level, these genes were also down-regulated, which was largely consistent with the results of bioinformatics analysis (Figure 9D).

Figure 9 Construction of the LPS-induced HUVECs and verification of gene expression. The expression of IL-6 (A) and TNF-α (B). (C and D) The expression of key PANoptosis-related sepsis predictor genes in HUVECs induced by LPS. ns: P ≥ 0.05; *P < 0.05; **P ≤ 0.01; ***p ≤ 0.001; ****P ≤ 0.0001.

Construction of a Nomogram Model and Survival Analysis

To evaluate the predictability of five key PANoptosis-related sepsis predictor genes, we further constructed a nomogram for risk assessment in sepsis samples (Figure 10A). The calibration curve indicated that the error between the actual and predicted risks was relatively small (Figure 10B). As demonstrated by the DCA, our nomogram has excellent value for clinical applications (Figure 10C). Correlation analysis between these genes and the results of the immune infiltration showed that all of these genes had a substantial negative correlation with CD4 naïve T cells and a strong positive correlation with activated dendritic cells and memory B cells (Supplementary Figure 1). Finally, we performed survival analysis on these genes, and the results showed that the Overall Survival of sepsis samples in the low-expression groups was almost all significantly higher compared with those in the high-expression groups (Figure 11).

Figure 10 Nomogram model of five key PANoptosis-related sepsis predictor genes. (A) Construction of a nomogram model. Calibration curve (B) and DCA (C) to evaluate the predictive performance.

Figure 11 Survival analysis of patients with different expression levels of (A) DDX60, (B) TAP2, (C) UBE2L6, (D) IFIT3, and (E) OAS2 in sepsis.

Discussion

Although several biomarkers of sepsis have been identified over the past decades, the early diagnosis of sepsis remains challenging due to the complexity of the intricate underlying pathophysiology and nonspecific early clinical manifestations. Therefore, it is essential to identify more suitable biomarkers to guide individualized sepsis treatment. Studies have shown that cell death does not proceed independently under pathological conditions but is often involved in extensive crosstalk.^23,24 PANoptosis is one such new concept emphasizing the interconnectedness among pyroptosis, apoptosis, and necroptosis, and that each of these cell deaths above individually shows a strong association with sepsis.^25–27 Nevertheless, the specific role of PANoptosis in sepsis remains inadequately investigated. Consequently, our research sought to investigate the specific function of PANoptosis in sepsis and utilize PANoptosis-related gene signatures to construct a gene prediction model.

Our present study firstly performed an in-depth analysis of the expression profiles of PRGs in sepsis. We observed that the expressions of PRGs were generally dysregulated in sepsis patients compared with normal samples, indicating a vital role of PRGs in the pathologic process of sepsis. Furthermore, immune infiltration analysis reveals that sepsis patients exhibited higher levels of infiltration of CD4+ T cells, monocytes, macrophages, and neutrophils, which aligns with findings from prior studies.^28–30 Additionally, we identified two distinct PANoptosis-related clusters in sepsis based on the expression of PRGs. DEGs related to these clusters were primarily enriched in pathways related to viral immune response, with notable activity in the NOD-like receptor signaling pathway and RIG-like receptors signaling pathway. Notably, NOD-like receptors (NLRs) and RIG-like receptors (RLRs) are essential members of the pattern recognition receptor (PRR) family, responsible for detecting invading microbial pathogens.³¹ Consistent with previously published studies on sepsis,^32–34 PANoptosis may indeed participate in the progression of sepsis through influencing innate immune pathways related to PRRs.

There has been a growing tendency to apply machine learning models for the prediction of sepsis,³⁵ which highlights that machine learning models take into account the relationships between clinical medicine and biological information, ultimately improving precision in sepsis diagnosis. Additionally, they provide insights into physiological pathways and potential therapeutic targets. In our current study, we performed a comparison of the forecast accuracy of four well-established machine learning models based on PANoptosis -specific DEGs. Then, five key PANoptosis-related sepsis predictor genes (DDX60, TAP2, UBE2L6, IFIT3, and OAS2) were finally selected to construct a 5-gene-based SVM model for further analysis, which exhibited satisfactory performance in anticipating the development of sepsis.

DDX60, belonging to the RNA helicase DEAD box (DDX) family, is an interferon (IFN)-inducible RNA helicase that has a useful part in activating RIG-I-like receptor-mediated signaling in the innate immune response.³⁶ Previous studies have shown that inhibiting DDX60 can significantly increase the level of CD4+ T cells, alleviate oxidative stress, and prolong the survival of septic mice.^37,38 The Transporter Associated with Antigen Processing 2 (TAP2) gene is a major histocompatibility complex (MHC) gene situated between HLA-DQ and HLA-DP.^39,40 TAP2 has been identified as a reliable marker for patients with a Mars1 endotype, consistently associated with 28-day mortality in sepsis cohorts,¹³ suggesting that TAP2 might be beneficial for diagnosing sepsis. Ubiquitin-conjugating enzyme E2 L6 (UBE2L6) is a specific enzyme that catalyzes ISGylation by conjugating interferon-stimulated gene 15 (ISG15) to target proteins⁴¹ and has a vital part in numerous pathophysiological processes, including responses to pathogen infections⁴² and immune responses.⁴³ Interferon‐inducible Protein with Tetrapeptide Repeats 3 (IFIT3) is closely correlated with a wide range of biological processes, which comprises apoptosis, differentiation, and cancer development,⁴⁴ and related studies have demonstrated that its absence can alleviate the inflammatory responses in mice with sepsis and improve survival outcomes in those with septic shock.^45,46 2’-5’-Oligoadenylate Synthetase 2 (OAS2) is an antiviral interferon-stimulated gene (ISG) and a member of the OAS protein family.⁴⁷ It has been reported that microvesicles extracted from the plasma in sepsis patients increase the level of OAS2 compared to healthy controls.⁴⁸ Collectively, these studies suggest that DDX60, TAP2, UBE2L6, IFIT3, and OAS2 may have significant implications in sepsis, leading to organ dysfunction, and our study further suggests that they may be associated with the regulation of PANoptosis in sepsis. Consequently, they hold promise as potential targets for therapeutic interventions in sepsis.

The high accuracy of five key PANoptosis-related sepsis predictor genes was validated on an external validation dataset (AUC = 0.989) in predicting sepsis. The mRNA and protein levels of nearly all genes, except for TAP2, were differentially expressed in sepsis and control samples, aligning with the findings from our previous bioinformatics analysis.

We also constructed a nomogram model for diagnosing sepsis subtypes using DDX60, TAP2, UBE2L6, IFIT3 and OAS2, which demonstrated significant predictive efficacy. Additionally, our observations indicated different expression of these genes was significantly associated with the prognosis of sepsis. Collectively, these findings suggest that the model could serve as a robust indicator for assessing the pathological outcomes and prognosis of sepsis patients.

Nevertheless, it’s important to acknowledge some limitations of our study. Firstly, we relied only on comprehensive bioinformatics analyses, which have not been clinically validated and may have some bias due to datasets from a single region. Therefore, further large clinical cohort studies are needed to confirm the predictive value of five key PANoptosis-related genes in sepsis. Secondly, additional experimental and clinical studies are needed to ascertain the precise mechanisms, by which these genes are involved in PANoptosis and inflammation in sepsis, as well as uncover the function of these genes in sepsis and explore potential early diagnostic biomarkers.

Conclusion

In summary, our study revealed dysregulations in PRGs and their associations with infiltrating immune cells, patient subtypes, and the optimal machine learning model for sepsis. Results from our bioinformatics analysis enhance our understanding of sepsis diagnosis and offer promising targets for its therapeutic interventions.