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/).

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.

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).

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.

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 .

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.

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).

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.

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.

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.

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 ).

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.

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.

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 ).

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 .

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.

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 ).

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).

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.

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 ).

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).