Construction and validation of a robust prognostic model based on immune features in sepsis

Purpose: Sepsis, with life-threatening organ failure, is caused by the uncontrolled host response to infection. Immune response plays an important role in the pathophysiology of sepsis. Immune-related genes (IRGs) are promising novel biomarkers that have been used to construct the diagnostic and prognostic model. However, an IRG prognostic model used to predict the 28-day mortality in sepsis was still limited. Therefore, the study aimed to develop a prognostic model based on IRGs to identify patients with high risk and predict the 28-day mortality in sepsis. Then, we further explore the circulating immune cell and immunosuppression state in sepsis.

Materials and methods: The differentially expressed genes (DEGs), differentially expressed immune-related genes (DEIRGs), and differentially expressed transcription factors (DETFs) were obtained from the GEO, ImmPort, and Cistrome databases. Then, the TFs-DEIRGs regulatory network and prognostic prediction model were constructed by Cox regression analysis and Pearson correlation analysis. The external datasets also validated the reliability of the prognostic model. Based on the prognostic DEIRGs, we developed a nomogram and conducted an independent prognosis analysis to explore the relationship between DEIRGs in the prognostic model and clinical features in sepsis. Besides, we further evaluate the circulating immune cells state in sepsis.

Results: A total of seven datasets were included in our study. Among them, GSE65682 was identified as a discovery cohort. The results of GSEA showed that there is a significant correlation between sepsis and immune response. Then, based on a P value <0.01, 69 prognostic DEIRGs were obtained and the potential molecular mechanisms of DEIRGs were also clarified. According to multivariate Cox regression analysis, 22 DEIRGs were further identified to construct the prognostic model and identify patients with high risk. The Kaplan–Meier survival analysis showed that high-risk groups have higher 28-day mortality than low-risk groups (P=1.105e-13). The AUC value was 0.879 which symbolized that the prognostic model had a better accuracy to predict the 28-day mortality. The external datasets also prove that the prognostic model had an excellent prediction value. Furthermore, the results of correlation analysis showed that patients with Mars1 might have higher risk scores than Mars2-4 (P=0.002). According to the previous study, Mars1 endotype was characterized by immunoparalysis. Thus, the sepsis patients in high-risk groups might exist the immunosuppression. Between the high-risk and low-risk groups, circulating immune cells types were significantly different, and risk score was significantly negatively correlated with naive CD4+ T cells (P=0.019), activated NK cells (P=0.0045), monocytes (P=0.0134), and M1 macrophages (P=0.0002).

Conclusions: Our study provides a robust prognostic model based on 22 DEIRGs which can predict 28-day mortality and immunosuppression status in sepsis. The higher risk score was positively associated with 28-day mortality and the development of immunosuppression. IRGs are a promising biomarker that might facilitate personalized treatments for sepsis.

Introduction

Sepsis is a complex disorder that develops as a severe systemic inflammatory response to infection, and is associated with high mortality (1). According to the US report, there were 48.9 (38.9-62.9) million incident cases of sepsis in the world annually and 11.0 (10·1-12·0) million patients died with sepsis (2). Sepsis was recognized as the most expensive burden and threat to human health. Increased mortality was associated with delay in initiating early treatments. The previous study estimated that the survival rate decreases by roughly 10% every hour that appropriate antimicrobial medication is delayed, emphasizing the urgent need for early identification and precise treatments to improve clinical outcomes (1, 3, 4). In 2017, World Health Organization (WHO) also declared that the improvement of sepsis early prevention, early recognition, and treatment is a global health priority (5). Therefore, the identification of septic patients at high risk may help clinicians to screen and identify individuals who are most likely to have poor prognosis, or to detect immunosuppressed states which could benefit from targeted immunostimulating therapies, and eventually improve patients’ prognosis.

Sepsis is an uncontrolled inflammatory response to invasive infection which can disturb homeostasis. After infection, the immune response can eliminate the pathogens but sometimes the host will release damage-associated molecular patterns (DAMPs) to damage organs. However, in late sepsis, sepsis patients have immune suppression which is characterized the lymphocyte exhaustion and the reprogramming of antigen−presenting cells (6, 7). In face of the complex pathophysiology of sepsis and its often challenging clinical evaluation, promising diagnostic biomarkers in sepsis are emerging with the application of blood genomics. Scicluna et al. (8) established endotypes for patients with sepsis through genome-wide blood gene expression profiles. The study provided a method to classify sepsis patients into four different endotypes and the detection of sepsis endotypes may assist in the precise treatments. Besides, increasing studies have identified novel immune biomarkers for early diagnosis and guide immunotherapies in oncology research. The immune related-genes (IRGs) model had been successfully applied in oncology to identify patients at high risk and estimate overall survival (9, 10). However, a robust IRGs model to identify high-risk patients and predict prognosis for adult patients with all-cause sepsis is still lacking.

The primary objective of this study was to construct an IRGs model to predict the prognosis of adult patients with all-cause sepsis. To achieve this aim, we obtain the differentially expressed IRGs (DEIRGs) that we can establish the Cox prediction model based on DEIRGs to predict the patients at high risk and the prognosis for patients with sepsis. Then, we construct a regulatory network between differentially expressed transcription factors (DETFs) and DEIRGs to explore the underlying molecular mechanisms. Besides, we further analyze the immune microenvironment in sepsis patients. Finally, we tested the robustness of the predictive model across the other datasets, and we provided a quantitative tool for predicting the individual probability of death.

Materials and methods

Datasets selection, data acquisition, and processing

A workflow is shown in Figure 1. The Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) and ArrayExpress (https://www.ebi.ac.uk/arrayexpress/) databases were comprehensively searched from inception to April 2022 to obtain the relevant datasets. The inclusion criteria of datasets were: (1) diagnosis of patients with sepsis; (2) sample size more than 50; (3) age ≥18 years; (4) the endpoints included 28-day mortality; (5) the patient’s specimens were collected before 24h on ICU admission and anti-inflammation treatments. Therefore, the 7 datasets were included in the study. Among them, GSE65682 was recognized as a training set because it was a large cohort study, and the other datasets were retained for model validation. The details of these datasets are shown in Table 1. To identify IRGs, we downloaded 2,483 IRGs from the Immunology Database and Analysis Portal (ImmPort) (http://www.immport.org/). Moreover, to construct the regulatory network, we obtained the transcription factors (TFs) from Cistrome Project (http://www.cistrome.org/).

FIGURE 1

Figure 1 Flowchart of data analysis and validation.

TABLE 1

Table 1 Basic information of the datasets included in this study.

Differential expression analysis in sepsis

All the genes in GSE65682 were differentially analyzed by using limma R packages (http://www.bioconductor.org/packages/release/bioc/html/limma.html) (11). The parameter for DEGs screened was│Log2Foldchange│≥0.5 and P-value < 0.05. The Volcano plots were drawn by ‘ggplot2’ R package. Then, the IRGs that were overlapping with DEGs were identified as DEIRGs. Similarly, DETFs were obtained by matching TFs with DEGs.

DEGs and gene set enrichment analysis

GSEA was used to assess related pathways and molecular mechanisms in sepsis. We performed the GSEA by the R package ‘clusterProfiler’. Normalized enrichment score (NES) and false discovery rate (FDR) were used to quantify enrichment magnitude and statistical significance, respectively (12, 13).

Identification of prognostic DEIRGs and construction of the regulatory network

To identify the prognostic DEIRGs (P <0.01), R package ‘survival R’ was used to perform univariate Cox regression analysis. Then, it is important to further explore the mechanisms of TFs to regulate the prognostic DEIRGs. Thus, we further analyzed the coexpression relationship between TFs and prognostic DEIRGs by calculating Pearson’s correlation coefficient. A regulatory network was constructed based on the filter thresholds (P value <0.001 and |cor| > 0.5). The network was visualized by using Cytoscape software.

Construction of the prognostic prediction model in sepsis and development of nomogram

Based on the univariate Cox regression analysis, prognostic DEIRGs were recognized as the biomarkers for multivariate Cox regression analysis. According to the median risk score value, conducted between low-risk and high-risk groups by using ‘survival’ R package. To evaluate the sensitivity and specificity of the prediction model, the receiver operating characteristics (ROC) curve was calculated using the ‘survivalROC’ package. The area under the ROC curve (AUC) was used to evaluate the prognostic model: 0.5-0.7 (moderate), 0.7-0.8 (better), and >0.9 (excellent).

To provide a quantitative tool for predicting probability of 28-day mortality in septic patients, we construct a nomogram according to the DEIRGs in the prognostic model and clinic features. The patients' clinic features are shown in Supplementary Table 1.

Validation in multiple external datasets

To evaluate the predictive performance of the prognostic model, 6 datasets (GSE63062, GSE95233, GSE106878, E-MTAB-4451, E-MTAB-5273, and E-MTAB-5274) were included according to the inclusion criteria. The prognostic model was used to predict the 28-day mortality of external datasets. Furthermore, the ROC curve was generated to determine sensitivity and specificity in the prognostic model.

Gene ontology and pathway enrichment analysis for DEIRGs in the prognostic model

To explore the mechanisms and functions of DEIRGs in the prognostic model, we performed Gene Ontology (GO) and Pathway Enrichment Analysis (KEGG) through the DAVID database (https://david.ncifcrf.gov/). Upon GO analysis and KEGG analysis, a P value <0.05 was recognized as statistical significance. The results of GO analysis were classified into three functional groups: biological process (BP), molecular function (MF), and cellular component (CC).

Correlation analysis between clinical features and DEIRGs in prognostic model

The correlation between risk score, gene expression value, and clinical features (age, gender, diabetes, ICU acquired infection (ICUA), and endotype class) were analyzed by using the ‘beeswarm’ R package. A P value <0.05 indicated statistical significance.

Exploration of circulating immune cells between low-risk and high-risk groups

CIBERSORTx (https://cibersort.stanford.edu/), an online analytical tool based on a kind of deconvolution algorithm iterated 1000 times, was available to provide an estimation of the abundances of member cell types in a mixed cell population by using gene expression data (14). Then, the content of 22 types of circulating immune cells in each sample was visualized by a vertical stack bar. Furthermore, the difference analysis of immune cells between low-risk and high-risk groups was shown by drawing barplot diagrams. Additionally, we explored the correlation between immune cells and risk score by Spearman correlation analyses. A P value <0.05 was considered statistical significance.

Statistical analysis

All statistical analyses were performed using R software and Grapad prism 9.0. The ‘limma R’ package was used to conduct differential expression analysis. The R package ‘clusterProfiler’ was adopted for assessing related pathways and molecular mechanisms in sepsis. The prognostic prediction model was constructed by univariate and multivariate Cox regression analysis. Besides, the ‘survival’, ‘survival ROC’, and ‘risk Plot’ R packages were applied to evaluate the survival difference between the high-risk and low-risk groups and assess the sensitivity and specificity in the prognostic model. Then, the ‘beeswarm’ R package was used to explore the correlation between clinical features and DEIRGs in the prognostic model. P value < 0.05 was considered statistically significant.

Results

DEGs, DEIRGs, DETFs and GSEA analysis

After the differential expression analysis of GSE65682, we obtained 3,648 DEGs (FDR <0.05, │Log2FC│≥0.5) (Figure 2A; Supplementary Table 2). To identify the DEIRGs, we downloaded all 2,483 immune genes from the ImmPORT database. Then, we matched IRGs with DEGs and obtained 278 DEIRGs (Figures 2B, D; Supplementary Table 3). Similarly, we searched and downloaded all 1,560 TFs from the Cistrome database. We matched TFs with DEGs and obtained 348 DETFs (Figures 2C, E; Supplementary Table 4).

FIGURE 2

Figure 2 Screening DEGs, DEIRGs and DETFs. (A) Volcano plot showing DEGs in GSE66890; (B) Venn diagram showed DEIRGs; (C) Venn diagram showed DETFs; (D) Volcano plot showing DEIRGs; (E) Volcano plot showing DETFs. Based on the |fold change|>0.5 and FDR<0.05, the red points represent upregulated genes and the green points represent downregulated genes. No significant differences are showed in black.

In order to explore the immune response in sepsis, we downloaded the KEGG gene sets and all GO gene sets from MsigDB. Then, the changes of these pathways and functions in gene sets were analyzed, namely Healthy vs Sepsis. According to our analysis, we found that immune response played an important role in the development of sepsis. In GO gene sets, we found that adaptive immune response was significantly upregulated in sepsis. However, cell activation involved in immune response, immune effector process, and myeloid leukocyte activation were upregulated in healthy and sepsis (Figures 3A, B). In KEGG gene sets, antigen processing and presentation, natural killer cell-mediated cytotoxicity, and primary immunodeficiency were most significantly increased (Figures 3C, D). Our results showed that there is a significant correlation between sepsis and immune response, and provide a theoretical basis for the construction of immune genes model to predict the prognosis of sepsis patients.

FIGURE 3

Figure 3 Exploring the difference of immune response between sepsis and healthy by using GSEA. (A) The enriched gene sets in GO collection; (B) The results of GO analysis from GSEA; (C) The enriched gene sets in KEGG collection; (D) The results of KEGG analysis from GSEA. The enrichment score of curve above 0 points indicates that the gene sets were activated in healthy. The curve below 0 points indicates that the gene sets were activated in sepsis. p.adjust, adjusted p-value; NES, normalized enrichment score.

Identification of prognostic DEIRGs and construction of regulatory network

Univariate Cox regression analysis was applied to screen and identify the prognostic genes in sepsis. According to P value <0.01, 69 prognostic DEIRGs were obtained (Figure 4A; Supplementary Table 5). Among them, 11 prognostic DEIRGs were high-risk and the others were low-risk. Then, to explore the molecular mechanisms between DETFs and prognostic DEIRGs, a regulatory network between DETFs and prognostic DEIRGs was constructed (Figure 4B; Supplementary Table 6). A total of 69 prognostic DEIRGs and 10 TFs were shown in the regulatory network (Figure 4B). As shown in the regulatory network, almost all expression level of high-risk DEIRGs (MPO, PTX3, DEFA4, CTSG, AZU1, ELANE, and RNASE3) were upregulated by CEBPE. Additionally, IL1R2 was upregulated by BCL11B, and FURIN was regulated by KLF1, TFDP1, and MX11. Besides, most low-risk DEIRGs had a positive relationship with BCL11B, MYC, POLB, STAT1, RUNX2, and KLF10. The other low-risk DEIRGs (HCK, IL17RA, ISG20L2, and ITGAL) were negatively regulated by KLF1, TFDP1, and MX11. The coefficient filter >0.5 and the P value <0.001 were set as the threshold to indicate statistical significance.

FIGURE 4

Figure 4 Prognostic DEIRGs and regulatory network between DETFs and prognostic DEIRGs. (A) Forest plot for prognostic DEIRGs in sepsis. Red and green dots were recognized as high-risk and low-risk, respectively; (B) Regulatory network between prognostic DEIRGs and DETFs. The red and green circles indicate high-risk DEIRGs and low-risk DEIRGs, respectively. The yellow triangles were applied to symbolize the DETFs. Moreover, the red and green lines were used to indicate a positive and negative correlation between prognostic DEIRGs and DETFs.

Construction of prognostic prediction model in sepsis

The 69 prognostic DEIRGs were obtained by univariate Cox regression analysis. Then, these prognostic DEIRGs were further incorporated into multivariate Cox regression analysis. Finally, 22 DEIRGs might serve to be the prognostic factors to independently predict the prognosis of sepsis patients (Table 2). Thus, the expression profiles of 22 DEIRGs were applied to construct the prognostic model to predict the 28-day mortality in sepsis patients. To obtain the survival risk score, the expression value and relative coefficients of 22 DEIRGs were used to calculate. The formulas was shown in Supplementary Table 7. Based on the median risk score value, 479 septic patients were classified into a high-risk group (n= 239) and a low-risk group (n=240) (Supplementary Table 7).

TABLE 2

Table 2 Multivariate Cox regression analyses of 22 IRGs of risk model in sepsis.

Then, Kaplan–Meier survival analysis was performed to analyze the 28-day mortality of high-risk groups (n= 239) and low-risk groups (n=240). As expected, the 28-day mortality of high-risk groups was significantly higher than low-risk group (Figure 5A; Supplementary Table 8). Furthermore, we drew an ROC curve to evaluate the sensitivity and specificity of the prognostic model. The results showed that the AUC value was 0.879 which symbolized that the prognostic model had a better accuracy to predict the 28-day mortality of high-risk and low-risk groups (Figure 5B). Additionally, the riskscope curve was constructed (Figure 5C) and the survival status of the two groups is shown in Figure 5D. The differential expression analysis of 22 DEIRGs are shown in Figure 5E.

FIGURE 5

Figure 5 Construction of prognostic model based on 22 DEGs. (A) Kaplan–Meier survival analysis of 28-day mortality between high-risk groups (red) and low-risk groups (blue). The color of each survival line indicated the 95% CI of probability of survival at each time point. (B) The ROC curve showed the AUC value of prognostic model. (C) The risk score analysis between high-risk and low-risk groups. (D) The survival status analysis between high-risk and low-risk groups. (E) The differentially expression analysis of 22 DEIRGs in prognostic model from 479 sepsis patients.

Validation of prognostic model by external datasets

To further evaluate the accuracy and reliability of the prognostic model, six datasets in line with the inclusion criteria were chosen to perform external validation. ROC analysis was performed to investigate the prognostic value of the prediction model. The AUC was 0.805 in E-MTAB-4451 (Figure 6A), AUC was 0.783 in E-MTAB-5273 (Figure 6B), AUC was 0.913 in E-MTAB-5274 (Figure 6C), AUC was 0.917 in GSE95233 (Figure 6D), AUC was 0.796 in GSE106878 (Figure 6E) and AUC was 0.915 in GSE63042 (Figure 6F), respectively. Therefore, the IRGs prognostic model had an excellent prediction value.

FIGURE 6

Figure 6 The prognostic efficacy of IRGs prognostic model. (A) The ROC curve of E-MTAB-4451 dataset. (B) The ROC curve of E-MTAB-5273 dataset. (C) The ROC curve of E-MTAB-5274 dataset. (D) The ROC curve of GSE95233 dataset. (E) The ROC curve of GSE106878 dataset. (F) The ROC curve of GSE63042 dataset.

Independent prognosis analysis and exploring the relationships between DEIRGs in prognostic model and clinical features in sepsis

The 22 DEIRGs in the prognostic model had a better predictive ability to investigate the 28-day mortality in sepsis. Then, we further conducted the univariate independent prognostic analysis and multivariate independent prognostic analysis to explore the correlation between clinical features and 28-day mortality in sepsis. The results of the univariate independent prognostic analysis showed that age (P = 0.019) and risk score (P <0.001) were related to the 28-day mortality, respectively (Supplementary Table 9; Figure 7A). The results of the multivariate independent prognostic analysis also showed that age (P = 0.005) and risk score (P <0.001) were the independent prognostic factors to predict the 28-day mortality in sepsis (Supplementary Table 10; Figure 7B).

FIGURE 7

Figure 7 The results of univariate independent prognostic analysis and multivariate independent prognostic analysis. (A) Univariate independent prognostic analysis. (B) Multivariate independent prognostic analysis. The red dots and green dots in the forest map indicated that the clinical feature was a high-risk factor and low-risk factor, respectively. ICUA, ICU acquired infection.

Then, the correlation between clinical features and DEIRGs in the prognostic model was further explored (Supplementary Material 11). Among the clinical features, the endotype class was classified into four classes, including Mars1, Mars2, Mars3, and Mars4 (8). According to the research, the Mars1 endotype was characterized by a pronounced decrease in the expression of genes corresponding to key innate and adaptive immune cell functions such as Toll-like receptor, nuclear factor κB (NFκB1) signaling, antigen presentation, and T-cell receptor signaling, which might be characterized by immune paralysis. The other endotypes (Mars2-4) were characterized by high expression of genes involved in pro-inflammatory (eg, NF-κB signaling) and innate (eg, interferon signaling) immune reactions, which are characterized as pro-inflammatory and innate immune response. As shown in Figure 8, the expression level of CD1D was higher in ICU-acquired infection (ICUA). Besides, the expression levels of ADRB2, CD1D, CD74, FYN, GNLY, IL16, IL17RA, PLXNC1, PSME1, TAP2, TNFRSF10B and TNFSF12 in Mars2-4 were significantly higher than Mars1. The expression levels of DEFA4, ELANE, MPO, and TFRC were significantly lower in Mars1 compared to those in Mars2-4. Therefore, DEFA4, ELANE, MPO, and TFRC might be related to immune paralysis in sepsis. Additionally, patients with Mars1 might have higher risk scores than Mars2-4 (Figure 8R) which was consistent with the results of Scicluna et al. (8).

FIGURE 8

Figure 8 Relationships between clinical features and DEIRGs in prognostic model. (A) Different expression of CD1D between the ICUA/NICUA in sepsis. (B–R) Different expression of DEIRGs in prognostic model between Mars1 and Mars2-4 in sepsis. ICUA, ICU acquired infection; NICUA, No ICU acquired infection; Mars, molecular diagnosis and risk stratification of sepsis.

Development of nomogram to predict the 28-day mortality in sepsis

We constructed a nomogram to predict the 28-day mortality in sepsis according to clinical features and DEIRGs in the prognostic model (Figure 9). The value of each of the variables was given a score based on the points scale axis. The total score was calculated by adding each single score. Then, the total points were projected to the 28-day mortality probability scale axis to estimate the probability of death in sepsis.

FIGURE 9

Figure 9 A constructed nomogram for 28-day mortality prediction of a patients with sepsis. ICUA, ICU acquired infection; Mars, molecular diagnosis and risk stratification of sepsis.

Functional analysis for DEIRGs in prognostic model

To explore the functional changes for DEIRGs in the prognostic model, we performed the functional enrichment analysis. The GO terms were divided into three functional groups, including biological process (BP), cell component (CC), and molecular function (MF). The top 10 significant enrichment results are shown in Figure 10. In BP groups, DEIRGs were mainly enriched in antigen processing and presentation, positive regulation of cytokine production and positive regulation of leukocyte cell−cell adhesion (Figure 10A). In CC groups, DEIRGs were mainly involved in MHC class II protein complex, MHC protein complex and phagocytic vesicle (Figure 10B). In MF groups, DEIRGs mainly enriched in cytokine binding, cytokine receptor activity and immune receptor activity (Figure 10C). As for the KEGG analysis, DEIRGs were mainly involved in antigen processing and presentation, cytokine−cytokine receptor interaction and hematopoietic cell lineage (Figure 10D).

FIGURE 10

Figure 10 The functional enrichment analysis for DEIRGs in prognostic model. (A) Biological process. (B) Cell component. (C) Molecular function. (D) KEGG pathway enrichment analysis.

Correlation analysis between DEIRGs and circulating immune cells

Numerous studies had demonstrated that circulating immune cells levels were associated with the prognosis of patients (15, 16). Therefore, we wanted to explore the different status of circulating immune cells between low-risk and high-risk groups. As shown in Supplementary Figure 1, the status of immune cells was significantly different in low-risk groups compared to the high-risk groups. Then, we further analyzed the composition of immune cells between low-risk and high-risk groups. The results of CIBERSORTx demonstrated that compared to the high-risk groups, CD8+ T cells (P=0.0135), resting (P=0.0005), and activated NK cells (P<0.0001), monocytes (P<0.0001), and M1 macrophages (P<0.0001) were more abundant in low-risk groups, while naive CD4+ T cells (P=0.0257), follicular helper T cells (P=0.0489) and activated dendritic cells (P<0.0001) were significantly enriched in high-risk groups (Figure 11A). Besides, we also analyzed the correlation of risk score and 22 immune cell types via Spearman correlation analyses. The results showed that risk scores were significantly positively correlated with follicular helper T cells (P=0.0437), gamma delta T cells (P=0.0004), resting NK cells (P=0.0259), activated dendritic cells (P<0.0001), and activated mast cells (P<0.0001), whereas were significantly negatively correlated with naive CD4+ T cells (P=0.019), activated NK cells (P=0.0045), monocytes (P=0.0134) and M1 macrophages (P=0.0002).

FIGURE 11

Figure 11 Comparison and correlation of circulating immune cells between low-risk and high-risk groups. (A) Comparison of circulating immune cells between low-risk and high-risk groups via CIBERSORTx. (B–J) Correlation between risk scores and circulating immune cells via Spearman correlation analysis. *p < 0.05.

Discussion

Sepsis, with high heterogeneity, is characterized by aberrant immune responses, including hyperinflammation and immune suppression (17). Increasingly, studies have pointed out that IRGs are promising novel biomarkers that may have important predictive and prognostic value (8, 18, 19). Thus, our research demonstrates that immune response played an important role in the development of sepsis. Then, the Cox prediction model obtained the 22 DEIRGs to classify the patients into low-risk and high-risk groups and construct the prognostic model. The regulatory network between TFs and prognostic DEIRGs was constructed to reveal the potential novel molecular mechanisms in sepsis. In this study, the prognostic model had a better accuracy to predict the 28-day mortality in sepsis. The external datasets also validated that the prognostic model had an excellent prediction value. Besides, we further developed a nomogram to provide a tool for predicting the probability of 28-day mortality in sepsis. We further explore the functional changes via functional enrichment analysis. Finally, the circulating immune cells were evaluated by CIBERSORTx.

As we know, the biomarkers to diagnose and predict the prognosis of sepsis were lacking due to the complex pathogenesis and high heterogeneity in sepsis. The unbalanced immune response of sepsis was initially activated to release tremendous damage-associated molecular patterns (DAMPs) such as cytokines. The cytokine storm will further lead to organ damage and even death (20). However, longitudinal analyses of immune response showed that patients developed persistent inflammation and immunosuppression in the late stage of sepsis (21). Therefore, the aberrant immune responses during sepsis might reflect the disease progression. The results of GSEA in this study also showed that immune responses were significantly related to the development of sepsis (Figure 3).

IRGs, promising novel biomarkers, had been used to predict the prognosis in many diseases (22). Even in sepsis, Lu et al. demonstrated the immune genes exhibited superior diagnostic and predictive efficacy in mortality than clinical characteristics (18). However, the molecular mechanisms of prognostic DEIRGs in sepsis were still unclear. In our study, we have identified 11 prognostic DEIRGs with high risk and 58 prognostic DEIRGs with low risk via univariate Cox regression analysis (Figure 4A). TFs, as an enhancer or promoter, could regulate the genes’ expression by binding to a particular DNA region. The regulatory network was constructed between TFs and prognostic DEIRGs (Figure 4B). We found key TFs (CEBPE, BCL11B, MYC, POLB, and STAT1) which had the most downstream DEIRGs and relatively high correlation coefficients. CEBPE has been demonstrated to be involved in the generation and proliferation of neutrophils (23). Besides, CEBPE was an important target to promote the innate immune system (e.g. neutrophil) to kill the bacteria (24). In this study, almost high-risk DEIRGs (MPO, PTX3, DEFA4, CTSG, AZU1, ELANE, and RNASE3) were upregulated by CEBPE which indicated that these high-risk DEIRGs might promote the inflammation and innate immune responses. Additionally, the BCL11B gene was essential for T cell and NK cell development and function (25). MYC, POLB, and STAT1 have also been described as having a strong relationship to the functions of the immune system and the clearance of pathogens (26–28). Therefore, most low-risk DEIRGs had a positive relationship with BCL11B, MYC, POLB and STAT1 might have an important role in regulating immune responses and defending against pathogens.

To construct the prognostic model, we further conducted the multivariate Cox regression analysis to identify the prognostic DEIRGs (Table 2). Then, the prognostic model could predict the 28-day mortality in sepsis with better accuracy (Figures 5B, 6). Among them, we obtained 4 high-risk genes (ELANE, IL1R2, MPO, and DEFA4) and 18 low-risk genes. ELANE and MPO had been demonstrated to involve in neutrophil protease activity. The expression levels of ELANE and MPO were correlated directly with organ failure and mortality which was in line with our results (29, 30). Besides, IL1R2, a decoy receptor for IL-1, has been implicated in sepsis (31). Previous studies have proven that IL1R2 was a biomarker to distinguish septic shock from non-septic shock postsurgical patients. The high expression of IL1R2 was significantly correlated to death in patients with postsurgical shock (32, 33). Interestingly, Liang et al. (34) pointed out that IL1R2 could distinguish gram-negative/gram-positive bacterial infection. The elevation of serum IL1R2 could be a biomarker to diagnose septic patients infected by gram-negative bacteria.

As we know, numerous studies have provided prognostic models/biomarkers for predicting overall survival in sepsis. However, these predictive factors were not applied to all sepsis patients due to the high heterogeneity. Thus, it is critical to stratify patients to guide treatments. A VANISH randomized trial categorized patients into SRS (sepsis response signatures) 1 and SRS2 according to transcriptomic profile. Patients with the immunocompetent SRS2 endotype might have significantly higher mortality when treated with corticosteroids than with placebo (35). Additionally, Scicluna et al. (8) classified patients with sepsis into four different endotypes (Mars1, Mars2, Mars3 and Mars4) upon ICU admission. According to the research, Mars1, with TAP2 transcripts denoting, was characterized by immune paralysis and poor prognosis, whereas Mars2-4 were characterized by high expression of pro-inflammatory genes. Our research also demonstrated that TAP2 was significantly downregulated in Mars1 compared to Mars2-4 (Figure 8N). TAP2 was a subunit of major histocompatibility complex class I (MHC-I) molecules involved in antigen processing (36). TAP2 has the potential to inhibit lipopolysaccharide-induced proinflammation by negative regulation of toll-like receptor-4 (TLR4) (37). Besides, our research also showed that most sepsis patients with high risk might have the Mars1 endotype which indicated the poor prognosis of patients with the Mars1 endotype (Figure 8R). This result was also in line with the previous study. Therefore, patients with sepsis in immunosuppression might be associated with an increased risk of mortality.

As we know, patients who survive early sepsis often develop a hypoinflammatory state and nosocomial infections which lead to high mortality (7, 17). Immune suppression in patients with sepsis is characterized by enhanced apoptosis of immune cells, T cell exhaustion, and reduced expression of activating cell surface molecules. Previous studies have proven that T cell exhaustion in immunosuppression was related to poor outcomes (38). The apoptosis of T cells (CD4+, CD8+, and Th17) will result in immunosuppression and is associated with higher mortality (39). Besides, nature killer (NK) cells could clear the pathogens and promote inflammation through the production of IFN-γ. However, NK cells will become tolerant and cytokine production of IFN-γ and TNF-α will be impaired in the late stage of sepsis. The proportion of NK cells in lymphocytes was negatively associated with 28-day mortality in septic patients (40, 41). Monocytes and macrophages are important components of the immune system that can remove pathogens and contribute to the immune response by antigen presentation. The M1 macrophages were characterized by the production of proinflammatory cytokines and antimicrobial activity. However, the polarization of M1 macrophages will be inhibited in immunosuppression. The M1 macrophage reprogramming will develop a pathological anti-inflammatory response to sepsis and increase the risk of immunosuppression (42, 43). In our research, CD8+ T cells (P=0.0135), resting (P=0.0005) and activated NK cells (P<0.0001), monocytes (P<0.0001), and M1 macrophages (P<0.0001) were more abundant in low-risk groups which indicated a hyperinflammatory state in low-risk groups (Figure 11A). Besides, the results of correlation analyses also showed that risk scores were significantly negatively correlated with naive CD4+ T cells (P=0.019), activated NK cells (P=0.0045), monocytes (P=0.0134), and M1 macrophages (P=0.0002) (Figure 11). In toto, the patients with high risk scores might be associated with immunosuppression. The risk score was positively associated with the development of immunosuppression in sepsis. The risk score might provide assistance for distinguishing sepsis patients with immunosuppression.

However, in spite of the remarkable results, there are several limitations that we could not ignore. First, a lot of publicly available sepsis datasets were excluded for lacking the mortality outcome. These datasets might concentrate on the differential diagnosis or other poor outcomes. The exclusion of these datasets might cause potential selection bias. Second, our prognostic model had a better performance in distinguishing patients with high risk, evaluating 28-day mortality in sepsis, and identifying sepsis patients with immunosuppression. However, it still needs large prospective cohorts to validate the performance before the prognostic model was applied to general use. Third, the datasets we included did not provide the details of comorbidities or other diseases. Therefore, we can’t exclude the impact of these factors on the prognostic model. Fourth, it may not accurately identify the immune cell types in sepsis according to bulk RNA-Seq data and the CIBERSORTx deconvolution algorithm. It still required further experiments (e.g. Flow Cytometry) to validate the results. Finally, the vivo and vitro experiments may help us identify the hub genes to predict the prognosis of sepsis and identify the patients with immunosuppression.

Conclusion

Our study demonstrated that immune response played an important role in the development of sepsis. IRGs, as promising novel biomarkers, were used to construct the TFs-DEIRGs regulatory network and prognostic prediction model, respectively. The TF-DEIRGs regulatory network has revealed the potential molecular mechanisms for DEIRGs in sepsis. The prognostic model, with great performance, could identify the patients with high risk and predict the 28-day mortality in patients with sepsis. Besides, the prognostic DEIRGs were also related to the immune cell circulating and immunosuppression state, which might promote individualized therapy for sepsis patients.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author contributions

YoZ and YH analysed and interpreted the data. YoZ, YuZ, YC and BL performed the bioinformatics analyses. XD, YX and LS performed data acquisition and figure preparations. YoZ and YH wrote the manuscript. XL and YL revised and edited this paper. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the National Natural Science Foundation of China (81870069, 81970071, 82070084, 82270085), the Natural Science Foundation of Guangdong Province (2020A1515011459, 2021A1515012565), the Science and Technology Program of Guangzhou (202102010366, 202201020444), the State Key Laboratory of Respiratory Disease Independent Program (SKLRD-Z-202108), Guangdong Marine Economy Development Special Project (GDNRC[2022]35).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.994295/full#supplementary-material

Supplementary Table 1 | The clinic features of GSE65682.

Supplementary Table 2 | The differential expression analysis of GSE65682.

Supplementary Table 3 | The differential expression analysis of IRGs.

Supplementary Table 4 | The differential expression analysis of TFs.

Supplementary Figure 1 | The status of circulating immune cells between low-risk and high-risk groups.

References

1. Singer M, Deutschman CS, Seymour CW, Shankar-Hari M, Annane D, Bauer M, et al. The third international consensus definitions for sepsis and septic shock (Sepsis-3). JAMA (2016) 315(8):801–10. doi: 10.1001/jama.2016.0287

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Rudd KE, Johnson SC, Agesa KM, Shackelford KA, Tsoi D, Kievlan DR, et al. Global, regional, and national sepsis incidence and mortality, 1990-2017: analysis for the global burden of disease study. Lancet (2020) 395(10219):200–11. doi: 10.1016/S0140-6736(19)32989-7

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Kumar A, Roberts D, Wood KE, Light B, Parrillo JE, Sharma S, et al. Duration of hypotension before initiation of effective antimicrobial therapy is the critical determinant of survival in human septic shock. Crit Care Med (2006) 34(6):1589–96. doi: 10.1097/01.CCM.0000217961.75225.E9

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Ferrer R, Martin-Loeches I, Phillips G, Osborn TM, Townsend S, Dellinger RP, et al. Empiric antibiotic treatment reduces mortality in severe sepsis and septic shock from the first hour: results from a guideline-based performance improvement program. Crit Care Med (2014) 42(8):1749–55. doi: 10.1097/CCM.0000000000000330

PubMed Abstract | CrossRef Full Text | Google Scholar

5. World Health Organization. World health assembly 70, resolution 70.7: improving the prevention, diagnosis and clinical management of sepsis (2017). Available at: http://apps.who.int/gb/ebwha/pdf_files/WHA70/A70_R7-en.pdf.

Google Scholar

6. van der Poll T, van de Veerdonk FL, Scicluna BP, Netea MG. The immunopathology of sepsis and potential therapeutic targets. Nat Rev Immunol (2017) 17(7):407–20. doi: 10.1038/nri.2017.36

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Cecconi M, Evans L, Levy M, Rhodes A. Sepsis and septic shock. Lancet (2018) 392(10141):75–87. doi: 10.1016/S0140-6736(18)30696-2

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Scicluna BP, van Vught LA, Zwinderman AH, Wiewel MA, Davenport EE, Burnham KL, et al. Classification of patients with sepsis according to blood genomic endotype: a prospective cohort study. Lancet Respir Med (2017) 5(10):816–26. doi: 10.1016/S2213-2600(17)30294-1

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Li B, Cui Y, Diehn M, Li R. Development and validation of an individualized immune prognostic signature in early-stage nonsquamous non-small cell lung cancer. JAMA Oncol (2017) 3(11):1529–37. doi: 10.1001/jamaoncol.2017.1609

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Pages F, Mlecnik B, Marliot F, Bindea G, Ou FS, Bifulco C, et al. International validation of the consensus immunoscore for the classification of colon cancer: a prognostic and accuracy study. Lancet (2018) 391(10135):2128–39. doi: 10.1016/S0140-6736(18)30789-X

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res (2015) 43(7):e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Eraso-Pichot A, Braso-Vives M, Golbano A, Menacho C, Claro E, Galea E, et al. GSEA of mouse and human mitochondriomes reveals fatty acid oxidation in astrocytes. Glia (2018) 66(8):1724–35. doi: 10.1002/glia.23330

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U.S.A. (2005) 102(43):15545–50. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods (2015) 12(5):453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Yang S, Liu T, Nan H, Wang Y, Chen H, Zhang X, et al. Comprehensive analysis of prognostic immune-related genes in the tumor microenvironment of cutaneous melanoma. J Cell Physiol (2020) 235(2):1025–35. doi: 10.1002/jcp.29018

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Song Q, Shang J, Yang Z, Zhang L, Zhang C, Chen J, et al. Identification of an immune signature predicting prognosis risk of patients in lung adenocarcinoma. J Transl Med (2019) 17(1):70. doi: 10.1186/s12967-019-1824-4

PubMed Abstract | CrossRef Full Text | Google Scholar

17. van der Poll T, Shankar-Hari M, Wiersinga WJ. The immunology of sepsis. Immunity (2021) 54(11):2450–64. doi: 10.1016/j.immuni.2021.10.012

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Lu J, Chen R, Ou Y, Jiang Q, Wang L, Liu G, et al. Characterization of immune-related genes andimmune infiltration features for early diagnosis, prognosis and recognition of immunosuppression in sepsis. Int Immunopharmacol (2022) 107:108650. doi: 10.1016/j.intimp.2022.108650

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Liu S, Li Y, She F, Zhao X, Yao Y. Predictive value of immune cell counts and neutrophil-to-lymphocyte ratio for 28-day mortality in patients with sepsis caused by intra-abdominal infection. Burns Trauma (2021) 9:tkaa040. doi: 10.1093/burnst/tkaa040

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Wiersinga WJ, Leopold SJ, Cranendonk DR, van der Poll T. Host innate immune responses to sepsis. Virulence (2014) 5(1):36–44. doi: 10.4161/viru.25436

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Darden DB, Kelly LS, Fenner BP, Moldawer LL, Mohr AM, Efron PA. Dysregulated immunity and immunotherapy after sepsis. J Clin Med (2021) 10(8):1742. doi: 10.3390/jcm10081742

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Li J, Liu C, Chen Y, Gao C, Wang M, Ma X, et al. Tumor characterization in breast cancer identifies immune-relevant gene signatures associated with prognosis. Front Genet (2019) 10:1119. doi: 10.3389/fgene.2019.01119

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Evrard M, Kwok IWH, Chong SZ, Teng KWW, Becht E, Chen J, et al. Developmental analysis of bone marrow neutrophils reveals populations specialized in expansion, trafficking, and effector functions. Immunity (2018) 48(2):364–379 e368. doi: 10.1016/j.immuni.2018.02.002

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Kyme P, Thoennissen NH, Tseng CW, Thoennissen GB, Wolf AJ, Shimada K, et al. C/EBPepsilon mediates nicotinamide-enhanced clearance of staphylococcus aureus in mice. J Clin Invest (2012) 122(9):3316–29. doi: 10.1172/JCI62070

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Holmes TD, Pandey RV, Helm EY, Schlums H, Han H, Campbell TM, et al. The transcription factor Bcl11b promotes both canonical and adaptive NK cell differentiation. Sci Immunol (2021) 6(57):eabc9801. doi: 10.1126/sciimmunol.abc9801

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Xu H, Jiang Y, Xu X, Su X, Liu Y, Ma Y, et al. Inducible degradation of lncRNA Sros1 promotes IFN-gamma-mediated activation of innate immune responses by stabilizing Stat1 mRNA. Nat Immunol (2019) 20(12):1621–30. doi: 10.1038/s41590-019-0542-7

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Wang S, Liu R, Yu Q, Dong L, Bi Y, Liu G. Metabolic reprogramming of macrophages during infections and cancer. Cancer Lett (2019) 452:14–22. doi: 10.1016/j.canlet.2019.03.015

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Alshalchi SA, Anderson GG. Involvement of stress-related genes polB and PA14_46880 in biofilm formation of pseudomonas aeruginosa. Infect Immun (2014) 82(11):4746–57. doi: 10.1128/IAI.01915-14

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Guirgis FW, Leeuwenburgh C, Moldawer L, Ghita G, Black LP, Henson M, et al. Lipid and lipoprotein predictors of functional outcomes and long-term mortality after surgical sepsis. Ann Intensive Care (2021) 11(1):82. doi: 10.1186/s13613-021-00865-x

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Almansa R, Heredia-Rodriguez M, Gomez-Sanchez E, Andaluz-Ojeda D, Iglesias V, Rico L, et al. Transcriptomic correlates of organ failure extent in sepsis. J infection (2015) 70(5):445–56. doi: 10.1016/j.jinf.2014.12.010

CrossRef Full Text | Google Scholar

31. Peters VA, Joesting JJ, Freund GG. IL-1 receptor 2 (IL-1R2) and its role in immune regulation. Brain Behav Immun (2013) 32:1–8. doi: 10.1016/j.bbi.2012.11.006

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Martinez-Paz P, Aragon-Camino M, Gomez-Sanchez E, Lorenzo-Lopez M, Gomez-Pesquera E, Fadrique-Fuentes A, et al. Distinguishing septic shock from non-septic shock in postsurgical patients using gene expression. J infection (2021) 83(2):147–55. doi: 10.1016/j.jinf.2021.05.039

CrossRef Full Text | Google Scholar

33. Martinez-Paz P, Aragon-Camino M, Gomez-Sanchez E, Lorenzo-Lopez M, Gomez-Pesquera E, Lopez-Herrero R, et al. Gene expression patterns distinguish mortality risk in patients with postsurgical shock. J Clin Med (2020) 9(5):1276. doi: 10.3390/jcm9051276

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Lang Y, Jiang Y, Gao M, Wang W, Wang N, Wang K, et al. Interleukin-1 receptor 2: A new biomarker for sepsis diagnosis and gram-Negative/Gram-Positive bacterial differentiation. Shock (2017) 47(1):119–24. doi: 10.1097/SHK.0000000000000714

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Antcliffe DB, Burnham KL, Al-Beidh F, Santhakumaran S, Brett SJ, Hinds CJ, et al. Transcriptomic signatures in sepsis and a differential response to steroids. from the VANISH randomized trial. Am J Respir Crit Care Med (2019) 199(8):980–6. doi: 10.1164/rccm.201807-1419OC

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Blees A, Januliene D, Hofmann T, Koller N, Schmidt C, Trowitzsch S, et al. Structure of the human MHC-I peptide-loading complex. Nature (2017) 551(7681):525–8. doi: 10.1038/nature24627

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Park S, Shin HJ, Shah M, Cho HY, Anwar MA, Achek A, et al. TLR4/MD2 specific peptides stalled in vivo LPS-induced immune exacerbation. Biomaterials (2017) 126:49–60. doi: 10.1016/j.biomaterials.2017.02.023

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Le Tulzo Y, Pangault C, Gacouin A, Guilloux V, Tribut O, Amiot L, et al. Early circulating lymphocyte apoptosis in human septic shock is associated with poor outcome. Shock (2002) 18(6):487–94. doi: 10.1097/00024382-200212000-00001

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Wu HP, Chung K, Lin CY, Jiang BY, Chuang DY, Liu YC. Associations of T helper 1, 2, 17 and regulatory T lymphocytes with mortality in severe sepsis. Inflammation Res (2013) 62(8):751–63. doi: 10.1007/s00011-013-0630-3

CrossRef Full Text | Google Scholar

40. Feng T, Liao X, Yang X, Yang C, Lin F, Guo Y, et al. A shift toward inhibitory receptors and impaired effector functions on NK cells contribute to immunosuppression during sepsis. J Leukoc Biol (2020) 107(1):57–67. doi: 10.1002/JLB.4A0818-313RR

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Souza-Fonseca-Guimaraes F, Parlato M, Fitting C, Cavaillon JM, Adib-Conquy M. NK cell tolerance to TLR agonists mediated by regulatory T cells after polymicrobial sepsis. J Immunol (2012) 188(12):5850–8. doi: 10.4049/jimmunol.1103616

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Essandoh K, Li Y, Huo J, Fan GC. MiRNA-mediated macrophage polarization and its potential role in the regulation of inflammatory response. Shock (2016) 46(2):122–31. doi: 10.1097/SHK.0000000000000604

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Rackov G, Hernandez-Jimenez E, Shokri R, Carmona-Rodriguez L, Manes S, Alvarez-Mon M, et al. p21 mediates macrophage reprogramming through regulation of p50-p50 NF-kappaB and IFN-beta. J Clin Invest (2016) 126(8):3089–103. doi: 10.1172/JCI83404

PubMed Abstract | CrossRef Full Text | Google Scholar