Transcriptomic data from human blood samples were obtained from the Gene Expression Omnibus (GEO) and ArrayExpress databases. Datasets containing samples from minors were excluded. Additionally, datasets that did not report mortality outcomes, lacked complete expression data, or included duplicated entries were also excluded. The expression data were obtained from public repositories and normalized according to the platform specifications using the limma package (Ritchie et al., 2015). Probes without corresponding gene symbols were removed. In cases where multiple probe sets mapped to the same gene symbol, the average expression value was used. Batch effects were corrected using the ComBat algorithm from the sva R package. Each dataset was treated as a batch, and expression matrices were merged on the intersection of shared genes after collapsing duplicate probes with limma. ComBat was applied in parametric empirical Bayes mode with mean–variance adjustment and batch indicators only, without additional covariates. The corrected expression matrix was used for downstream analyses (Leek et al., 2012; Zhang et al., 2020a). Only genes common across all included datasets were retained for model development and validation.

In this study, NMF was employed for unsupervised clustering based on Leptin-related genes from the GeneCard database (Stelzer et al., 2016), with the aim of identifying molecular subtypes of sepsis. Gene expression data for key genes were extracted from the GeneCard database, and the data were transposed to meet the input format required for NMF analysis. NMF was performed using the NMF package, with ranks ranging from 2 to 10 to explore different potential clusters of sepsis patients. The Brunet method was applied for matrix factorization.

Single Sample Gene Set Enrichment Analysis (ssGSEA) (Subramanian et al., 2005) was used to assess immune cell infiltration in blood samples from septic patients. This method estimates the enrichment of immune cell types in individual samples based on predefined gene sets corresponding to various immune cell types. The analysis was performed using the GSVA package (Hänzelmann et al., 2013), which calculates the enrichment score for each immune cell type in each sample. The enrichment score represents the relative abundance of each immune cell type within the respective sample.

Using the WGCNA package (Langfelder and Horvath, 2008), potential Leptin-related gene modules associated with sepsis were identified. The data were preprocessed by filtering out genes with low variance, retaining only those with variance greater than the median. Missing values in the dataset were checked, and genes or samples with excessive missing data were removed to ensure the integrity of the analysis. Hierarchical clustering was performed on the samples to detect potential outliers, with samples outside a predefined cutoff being excluded. Gene modules were identified using hierarchical clustering based on the Topological Overlap Matrix (TOM). The optimal soft-thresholding power was determined, with a scale-free topology criterion set at 0.9. The TOM was then calculated, followed by hierarchical clustering of genes to construct a dendrogram, which was used for module detection. Modules were identified using the dynamic tree cutting method, with a minimum module size of 50 genes. The identified modules were correlated with clinical traits, such as survival data, using Pearson’s correlation. Significant modules were identified based on the correlation between gene module membership (MM) and gene trait significance (GS). Finally, the relationship between gene expression and clinical traits was assessed by visualizing MM and GS values through scatter plots. Genes most closely related to significant modules were retained for further analysis, including the identification of potential biomarkers associated with sepsis and mortality.

To explore the biological functions and pathways associated with the identified gene modules and genes, Gene Ontology (GO) (The Gene Ontology Consortium, 2019) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000) functional enrichment analyses were conducted. GO analysis classifies genes into three main categories: biological process, molecular function, and cellular component. KEGG pathway enrichment analysis was then performed to identify the biological pathways involved.

DEGs were identified by applying a fold change (FC) cutoff of 0.8 and a p-value threshold of 0.05. The limma package (Ritchie et al., 2015) was used to perform the differential expression analysis, with p-values adjusted for multiple testing using the Benjamini-Hochberg method.

In this study, multiple machine learning algorithms were employed to construct prognostic classification models for predicting 28-day mortality in sepsis. These included combinations such as Lasso Regression + Random Forest (RF), Generalized Linear Model Boosting (GLMB) + RF, Lasso Regression + Gradient Boosting Machine (GBM), RF, GBM, and combinations of Stepwise Logistic Regression (SLR) with RF, among others. These models were trained on internal datasets and validated externally. Variable selection was performed using methods such as Lasso Regression and Elastic Net (EN). The performance of the models was evaluated using several metrics, including Receiver Operating Characteristic (ROC) curves, confusion matrix (CM), and decision curve analyses, to assess predictive accuracy, classification performance, and clinical utility. ROC curves assessed model discrimination, the CM evaluated classification capability, and DCA measured the potential clinical benefit for predicting 28-day mortality and informing decision-making.

In this study, various machine learning algorithms were employed to construct prognostic models for sepsis, including combinations of Random Survival Forest (RSF), Lasso Regression, Stepwise Cox Regression (StepCox), Gradient Boosting Machine (GBM), CoxBoost, Elastic Net (Enet), Ridge Regression, and others. The models were trained and evaluated using external validation datasets. To assess model performance, the concordance index (c-index) was calculated to evaluate the discriminatory power of the models in predicting patient outcomes. Additionally, time-dependent receiver operating characteristic (ROC) curves were employed to assess the predictive accuracy of the models over time, particularly for 28-day mortality prediction. The area under the curve (AUC) of the time-dependent ROC was computed to measure the ability of the models to classify patients correctly at various time points, providing a more comprehensive evaluation of their performance.

Single-cell RNA sequencing data for this study were retrieved from the GEO database under accession number GSE167363 (Qiu et al., 2021), encompassing transcriptomic profiles of peripheral blood mononuclear cells (PBMCs) derived from healthy individuals and patients with gram-negative sepsis (both survivors and non-survivors). Raw count matrices generated via the 10X Genomics platform were imported and processed using the Seurat (Butler et al., 2018; Slovin et al., 2021) framework in R. To retain high-confidence single-cell profiles, we applied stringent quality control filters: cells with fewer than 1,000 total RNA molecules, fewer than 200 or more than 10,000 detected genes, over 20% mitochondrial gene expression, or over 20% ribosomal gene content were excluded. This filtering strategy helped eliminate potential artifacts such as dead cells, doublets, or empty droplets. The resulting curated dataset was stored for subsequent integrative and functional analyses. Following rigorous quality assessment, gene expression values were normalized using the LogNormalize approach in the Seurat package, wherein each gene count was scaled relative to the total expression per cell and log-transformed after multiplication by a scale factor of 10,000. To capture transcriptional heterogeneity, highly variable genes (HVGs) were selected using the variance-stabilizing transformation (vst) method, retaining the top 2,000 genes. These features were subsequently standardized to zero mean and unit variance prior to dimensionality reduction. Principal component analysis (PCA) was carried out on the HVGs, and inter-sample variation was mitigated using Harmony (Korsunsky et al., 2019), leveraging sample identity as the integration variable. The corrected principal components were then used to generate a t-distributed stochastic neighbor embedding (t-SNE) projection for visualization. Clustering was performed using a graph-based shared nearest neighbor (SNN) approach, with resolution tuning informed by hierarchical visualization through clustree. Final clustering was executed at a resolution of 1.2. Cell type annotation combined automated classification via SingleR (Zhao et al., 2020)—using reference transcriptomes of human immune cells—with manual refinement based on canonical marker gene expression patterns.

A total of nine peripheral blood samples were obtained from Sichuan Provincial People’s Hospital, comprising three healthy controls, three sepsis patients who survived beyond 28 days, and three sepsis patients who died within 28 days of hospital admission. All samples were collected at the time of admission, prior to the initiation of any therapeutic intervention. Peripheral blood mononuclear cells (PBMCs) were isolated from freshly drawn blood using Ficoll-Paque density gradient centrifugation. The study protocol was approved by the Medical Ethics Committee of Sichuan Provincial People’s Hospital (Approval No. 2023-581), and written informed consent was obtained from all participants or their legal guardians.

Total RNA was extracted from PBMCs using Trizol reagent (Invitrogen, Carlsbad, CA, USA). cDNA synthesis was performed using the ReverTra Ace qPCR RT Kit (TOYOBO, Osaka, Japan) according to the manufacturer’s instructions. For real-time quantitative PCR detecting system (RQ-PCR) analysis, TransStart Tip Green qPCR SuperMix (TransGen Biotech, Beijing, China) was used. Gene expression was normalized to GAPDH, and relative expression was calculated using the ΔΔCt method. The primers are listed in Supplementary Table 1 .

PBMCs were lysed in RIPA buffer (Beyotime Biotechnology, Shanghai, China) containing 1 mM PMSF. Protein concentrations were measured using a BCA assay kit (Beyotime Biotechnology, Shanghai, China). Equal protein aliquots (30–50 μg) were separated by 8–12% SDS-PAGE and transferred to PVDF membranes (Servicebio, Wuhan, China). The membranes were blocked with 5% non-fat milk in TBST for 1 hour at room temperature, followed by overnight incubation at 4 °C with primary antibodies: PILRA (1:1,000), TFRC (1:1,000), and GAPDH (1:3,000) for loading control. After washing with TBST three times, membranes were incubated with HRP-conjugated goat anti-rabbit or anti-mouse secondary antibodies (1:3,000; Servicebio, Wuhan, China) for 1 hour at room temperature. Protein bands were visualized using ECL reagent (Servicebio) and analyzed with ImageJ software. GAPDH was used as a loading control.

The analysis in this study was performed using R software (version 4.3.1). Spearman’s rank correlation coefficient was used to evaluate correlations between variables. To assess survival differences among groups, Kaplan-Meier survival curves were generated and compared using the log-rank test. The optimal threshold for survival analysis was determined using the surv_cutpoint function, which identifies the cutoff point that maximizes the survival differences. For statistical comparisons, an independent t-test or Mann–Whitney U test was applied for two-group analyses, depending on data distribution, while one-way ANOVA or Kruskal–Wallis test was employed for comparisons among multiple groups. A P-value of less than 0.05 (two-sided) was considered statistically significant.

The basic information of the datasets used in this study is provided in Supplementary Table 2 . The code used in this study is provided in Supplementary Data Sheet S1. Genes associated with Leptin were retrieved from the GeneCards database ( Supplementary Table 3 ). Using these genes, NMF was performed on sequencing data from GSE65682. As shown in Figure 1A , the consensus matrix divided the samples into three distinct clusters. To determine the optimal number of clusters, an NMF rank survey was conducted ( Figure 1B ). Various evaluation metrics, including cophenetic, dispersion, evar, residuals, RSS (Residual Sum of Squares), silhouette, and sparseness, were assessed across factorization ranks from 2 to 10. The NMF rank survey revealed that three clusters provided the best fit according to multiple metrics, indicating the clearest separation between the clusters. PCA further supported the distinct separation of the three clusters ( Figure 1C ). Prognostic analysis showed significant differences among the clusters, with Cluster 1 exhibiting the best prognosis, while Cluster 3 had the worst outcome ( Figure 1D ). Furthermore, immune infiltration differences between the clusters were compared. With the exception of Gamma.delta. T.cells, significant statistical differences in immune cell infiltration were observed between the three clusters ( Figure 1E ). In the analysis of clinical characteristics in relation to the clusters, several clinical variables were evaluated for their distribution across the clusters ( Figures 1F-J ). Age did not differ significantly between clusters, with similar age distributions observed across the three groups. The gender distribution was balanced across the clusters, with no significant differences in the male-to-female ratio within each cluster. For thrombocytopenia, there were no notable differences across the three groups, with both low and normal platelet levels present in all clusters. The proportion of patients with ICU-acquired infections was similar across the clusters, with no significant variation. Additionally, the prevalence of diabetes mellitus was evenly distributed across all clusters, showing no significant differences.

Prior to performing WGCNA, the expression matrix derived from the GSE65682 dataset was confirmed to be of high quality, with no apparent outliers or missing values. A soft-thresholding power of β = 9 was selected, as it was the minimal value at which the scale-free topology model began to stabilize, as illustrated in Figure 2A . This threshold was subsequently used to generate the TOM and conduct preliminary module detection. Modules displaying similar gene expression profiles were further combined based on their eigengenes, yielding seven distinct gene modules ( Figure 2B ). Among them, the pink and salmon modules demonstrated particularly strong correlations, as shown in Figure 2C . Further analysis revealed significant associations between MM and GS within these modules ( Figures 2D, E ).

Functional enrichment analysis revealed distinct biological roles for genes in the pink and salmon modules ( Figures 2F, G ). Genes in the pink module were mainly involved in immune-related processes, including leukocyte migration, chemotaxis, and phagocytosis. They were enriched in components such as secretory granule membranes and focal adhesions, and showed functions related to GTPase regulation and kinase activity. KEGG pathways indicated involvement in phagocytosis, chemokine signaling, and immune cell migration. In contrast, genes in the salmon module were associated with cell cycle–related processes like chromosome segregation and nuclear division. These genes were enriched in vesicle lumen and lysosome components, with functions including DNA helicase activity and microtubule binding. KEGG enrichment highlighted pathways such as DNA replication, cell cycle, and mismatch repair, suggesting roles in cell proliferation and genomic stability.

To reduce technical variation among datasets, the ComBat algorithm was employed for batch effect correction. As shown in Figure 3A , prior to adjustment, samples from GSE54514, GSE65682, and GSE95233 clustered separately in the PCA plot, indicating substantial batch effects. After correction ( Figure 3B ), the datasets were well integrated, with samples from different cohorts intermixed, confirming effective removal of batch effects. Following integration, a total of 745 DEGs were identified based on predefined criteria ( Figure 3C ). These DEGs were intersected with gene sets from the pink and salmon modules derived from WGCNA, resulting in 70 overlapping hub genes ( Supplementary Table 4 ). Chromosomal mapping showed that these hub genes were distributed across most chromosomes, excluding chromosomes 4, 12, 18, and the Y chromosome ( Figures 3D, E ). Functional enrichment analysis revealed that these hub genes are primarily involved in immune-related biological processes such as bacterial defense, cell–cell adhesion, and response to Gram-negative bacteria ( Figure 3F ). GO analysis also highlighted their localization to secretory granule lumens and involvement in kinase regulation and lipopolysaccharide binding. KEGG pathway analysis ( Figure 3G ) showed enrichment in pathways including axon guidance, lysosome, efferocytosis, and folate biosynthesis, suggesting potential roles in immune modulation and cellular clearance mechanisms in sepsis.

A variety of machine learning algorithms were utilized to construct a prognostic classification model aimed at predicting 28-day mortality in sepsis. The GSE65682 dataset was designated as the training cohort, while E-MTAB-4451, E-MTAB-5273, E-MTAB-7581, and GSE63042 served as independent validation cohorts. In total, 113 models—including both individual and ensemble algorithm strategies—were assessed, and their respective AUC values are presented in Figure 4A . Among these, the Lasso combined with Random Forest (RF) approach achieved the highest mean AUC across datasets and was therefore selected to construct the final predictive model ( Supplementary Table 5 ) and compute the corresponding risk scores. The performance of the established risk score model was then evaluated in both the training and external validation sets using receiver operating characteristic (ROC) curves, confusion matrices, and clinical decision curve analysis (DCA) ( Figures 4B–F ). In the training set (GSE65682, Figure 4B ), the model exhibited strong discriminatory capacity with an AUC of 0.94, along with high sensitivity, specificity, and net clinical benefit. In external datasets including E-MTAB-5273 ( Figure 4C ), E-MTAB-7581 ( Figure 4D ), GSE63042 ( Figure 4E ), and E-MTAB-4451 ( Figure 4F ), the model retained stable and reliable performance, with AUC values ranging from 0.74 to 0.83. The confusion matrices demonstrated balanced classification accuracy, while the DCA curves suggested superior clinical utility of the model compared to conventional reference strategies across a wide range of threshold probabilities. Besides, ROC curve analysis in the GSE65682 cohort demonstrated that the proposed model outperformed traditional clinical variables in predictive performance ( Supplementary Figure 1A ).

Multiple machine learning algorithms were applied to develop a prognostic model for sepsis patients. The GSE65682 dataset was used as the training cohort, while GSE54514 and GSE95233 served as external validation sets. A total of 98 models, constructed using either single algorithms or algorithmic combinations, were evaluated based on their concordance index (C-index) across datasets ( Figure 4A ). Among them, the RSF + GBM combination achieved the highest average C-index and was selected to build the final prognostic model ( Supplementary Table 5 ) and compute individual risk scores ( Figure 5A ). As shown in Figures 5B, C , patients classified into the high-risk group had significantly worse survival outcomes compared to those in the low-risk group, with the exception of the GSE54514 dataset ( Figure 5D ). Furthermore, time-dependent ROC curves demonstrated that the risk score maintained strong prognostic accuracy across multiple time points. Finally, the model’s predictive performance was benchmarked against previously published prognostic models (Liang et al., 2022; Zheng et al., 2022; Chen et al., 2023; Lin et al., 2023; Liu et al., 2023; Chen et al., 2024) ( Figure 5E ). In the GSE65682 dataset, the proposed model outperformed all others, and in the remaining two cohorts, it consistently ranked among the top-performing models. Besides, C-index comparisons at 7, 14, and 28 days consistently demonstrated that the prognostic model achieved superior discriminative ability compared with conventional clinical variables ( Supplementary Figures 1B ). Moreover, both univariate and multivariate Cox regression analyses confirmed that the model-derived risk score was an independent prognostic factor for sepsis outcomes ( Supplementary Figures 1C, D ). Among them, five features—TFRC, PILRA, DEFA4, KRT23, and BEX1—are shared between the diagnostic and prognostic models.

After rigorous quality filtering to remove doublets, apoptotic cells, and empty droplets, high-confidence single-cell transcriptomic data were obtained ( Supplementary Figure 2A ). To correct sample-derived variability, the Harmony algorithm was applied, resulting in improved integration of cells across conditions ( Figure 6A ), with separate visualizations for the sepsis and control groups provided in Figure 6B . Clustering based on Harmony-corrected principal components (PC) (with PC 13 selected; Supplementary Figures 2B, C ) identified 23 transcriptionally distinct clusters ( Figure 6C ), which were annotated using the SingleR package and grouped into six major immune cell types: Monocyte, NK_cell, Neutrophils, Platelets, Erythroblast, and B_cell ( Supplementary Figure 2D ). Within the monocyte lineage, three subsets were identified ( Supplementary Figure 2E ): classical, intermediate, and non-classical (Wong et al., 2011; Thomas et al., 2017; Li et al., 2024). In total, eight immune cell populations were delineated, including Monocytes (subdivided into classical, intermediate, and non-classical subsets), NK cells, Neutrophils, Platelets, Erythroblasts, and B cells ( Figure 6D ). Among the five candidate genes (TFRC, PILRA, DEFA4, KRT23, and BEX1), TFRC and PILRA displayed significantly higher expression in monocyte subsets compared with other immune cell types, whereas DEFA4, KRT23, and BEX1 showed minimal or negligible expression across all populations. Feature plots in Figure 6E confirmed their cell-type-specific expression patterns, and both TFRC and PILRA were differentially expressed between the sepsis and control groups ( Figure 6F ). Notably, in the sepsis group, TFRC and PILRA also exhibited significant variation in expression across the eight immune cell types, whereas no such differences were observed for the remaining genes ( Figure 6G ).

Validation in PBMCs from healthy controls, sepsis survivors, and non-survivors confirmed that TFRC expression progressively increased, whereas PILRA decreased across these groups (qPCR, Figure 7A ; WB, Figures 7B, C ). Analysis of the GSE154401 dataset showed that leptin antibody treatment of T cells significantly elevated TFRC, while PILRA exhibited a non-significant reduction ( Figure 7D ). Consistently, in GSE57065 and GSE131761, TFRC was upregulated and PILRA was downregulated in sepsis compared with controls ( Figures 7E, F ).