The Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo) was used to access the GSE167363 (n = 12) and GSE175453 (n = 9) scRNA-seq sets. These two datasets (GSE167363 and GSE175453) utilized in this study can be found in Supplementary Tables 1 and 2 . For bulk RNA-seq data, we utilized the GSE57065 dataset (n = 107) as the training dataset and the GSE95233 dataset (n = 124) as the test dataset. These two datasets (GSE57065 and GSE95233) utilized in this study can be found in Supplementary Tables 3 and 4 . We selected OS-related genes (n = 807) from the Genecards database (https://www.genecards.org/). Genes are listed in Supplementary Table 5 .

In the processing of scRNA-seq data, we retained high-quality cells with less than 20% mitochondrial gene content and expressing more than 200 genes. Additionally, we prioritized genes with expression levels ranging from 200 to 7,000 and active in at least three cells. A total of 103,851 eligible cells were retained for further analysis. Subsequently, data integration was performed using the Seurat pipeline. The remaining cells were then scaled and normalized via a linear regression model employing the “Log-normalization” method. Additionally, the top 3,000 highly variable genes were identified using the “FindVariableFeatures” function. Following this, the dimensionality of the scRNA-seq data was reduced using principal component analysis (PCA). To eliminate batch effects among samples, soft k-means clustering was performed with the “Harmony” package. Cell clustering was subsequently carried out using the “FindClusters” function, with the resolution parameter adjusted to 0.6. The annotation of cell clusters was based on genes with high expression levels, genes displaying distinct expression patterns, and established canonical cell markers (Yang et al., 2024).

Five algorithms (including AUCell, Ucell, singscore, ssgsea, and AddModuleScore) were used to evaluate each cell’s OS activity at the single-cell level and determine the overall OS activity. The raw score matrix underwent sequential Z-score standardization and Min–Max normalization: initial standardization transformed features to zero-mean and unit-variance distributions, followed by a custom normalization function linearly scaling values to the [0,1] range applied column-wise to ensure feature-wise consistency. The processed matrix was converted into a data frame format, with a composite score derived from row-wise summation of normalized feature values. Based on the quartile method, cells with a score below the 25th percentile were categorized into the low-OS activity state group, those between the 25th and 75th percentiles were classified as those belonging to the “transition state”, and those above the 75th percentile were designated as the high-OS state group. The correlation analysis was used to explore the association between genes and OS activity. The “FindMarkers” function was used to perform differentially expressed gene (DEG) analysis in the upregulated OS group (avg_log2FC > 0.25, p _adj < 0.05). To assess whether the distribution of cell types differs significantly across the high-, medium-, and low-scoring groups, we constructed contingency tables for each cell type and performed chi-squared tests of independence. For cell types with small sample sizes, robustness was checked with Fisher’s exact test. The resulting p-values were adjusted for multiple testing using the false discovery rate (FDR) method.

The hdWGCNA was used to construct a co-expression network using the “hdWGCNA” package. We used the genes expressed in more than 5% of cells, and the “MetacellsByGroups” function was used to construct the metacell gene expression matrix. The soft power was determined by the “TestSoftPowers” function. The “ConstructNetwork” function is employed to build the co-expression network.

To identify optimal feature genes associated with overall survival, we performed an integrative analysis of the pre-screened gene set using five distinct machine learning approaches: (1) least absolute shrinkage and selection operator (LASSO) regression, (2) support vector machine–recursive feature elimination (SVM-RFE), (3) Boruta feature selection, (4) random forest (RF), and (5) gradient boosting machine (GBM). This multi-algorithm framework ensured robust identification of hub genes while mitigating model-specific biases. Feature genes were determined by intersecting the outputs of all five algorithms (LASSO, SVM-RFE, Boruta, RF, and GBM), ensuring robust selection. LASSO algorithms eliminated redundant predictors while maintaining discriminative power. SVM-RFE was employed to iteratively prune the feature set by removing the least informative features, thereby improving the model’s predictive performance. The Boruta algorithm assessed the significance of each feature by repeatedly sampling from the original dataset and constructing RFs. Meanwhile, the RF algorithm built multiple decision trees through random sampling and feature selection, and generated predictions via voting or averaging mechanisms. These approaches collectively clarified the correlations and interactions among features. GBM is a powerful ensemble technique that builds models iteratively to minimize a loss function. It is highly effective for a wide range of tasks but requires careful tuning and regularization to avoid overfitting.

We leveraged CellChat to analyze variations in cell–cell communication modules by integrating gene expression data. Following the standard CellChat pipeline, we employed the default CellChatDB as the ligand–receptor interaction database. Cell type-specific interactions were inferred by identifying overexpressed ligands or receptors within each cell cluster, and validated through the assessment of enhanced ligand–receptor interactions upon their upregulation. This approach ensured robust inference of communication patterns while maintaining alignment with established analytical frameworks. Additionally, the activity of gene regulatory networks was inferred using the R package “Scenic”, while for pseudotime analysis, three algorithms [“Monocyte” (Xu et al., 2024), slingshot (Street et al., 2018), and “CytoTRACE” (Wang et al., 2022)] were used for analysis.

The THP-1 human monocytic leukemia cell line (Type Culture Collection of Chinese Academy of Sciences, China) was cultured in RPMI-1640 medium (Gibco, USA) containing 10% fetal bovine serum (FBS) (Gibco, USA) with 2 mM L-glutamine (Solarbio, China), 1 mM sodium pyruvate (Solarbio, China), 10 mM HEPES (Gibco, USA), and penicillin/streptomycin (Solarbio, China) (50 U/mL each).

The differential protocol of macrophage was according to a previous study (Lai et al., 2023). Cells were maintained at 37°C in a 5% CO₂ humidified incubator. For lipopolysaccharide (LPS) (Selleck, USA) stimulation experiments, triplicate cultures (1×10⁶ cells/mL/well in 12-well plates) were treated with 10 ng/mL LPS for 24 h. Post-treatment cell pellets obtained by centrifugation (1,000 rpm, 5 min) were divided for parallel processing: RNA extraction using Trizol (TaKaRa Bio Inc., Japan) and protein isolation.

Gene silencing of LILRA5 was achieved using siRNA technology with Lipofectamine™ 2000 (Invitrogen, USA). Specifically, the siRNA duplexes targeting the sequence were commercially synthesized (Sangon Biotech, Shanghai). When cells reached 30%–50% confluency, transfection was performed with either LILRA5-targeting siRNA or non-targeting control (NC) using the Lipofectamine™ 2000 system. Following 6 h of transfection, cells were processed for downstream analyses.

RNA concentration was measured using the NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, MA, USA). For reverse transcription, the SuperScript First-Strand Synthesis System RT-PCR kit (TaKaRa Bio Inc., Japan) was employed to synthesize complementary DNA (cDNA) from RNA. Quantitative PCR was conducted using the GoTAq qPCR Master Mix (Promega, USA) on the Rotor-Gene Q PCR detection system (Qiagen, Germany). A 10-µL reaction mixture was subjected to thermal cycling, beginning with an initial denaturation at 95°C for 10 min. This was followed by 40 cycles of amplification, consisting of 95°C for 5 s, 60°C for 30 s, and 72°C for 30 s. Subsequently, gene expression levels were quantified using the2^–ΔΔCt method. Primer sequences are listed in Table 1 .

THP-1 cells were seeded at a density of 1×10⁶ cells/mL in triplicate wells of 12-well plates. Experimental groups were stimulated with LPS (10 ng/mL) for 30 min, while control groups were treated with phosphate-buffered saline (PBS) (Thermo Fisher, USA). Proteins were extracted from THP-1 cells using RIPA (Solarbio, China) and phenylmethanesulfonyl fluoride (PMSF) (Boster, USA). After polyacrylamide gel electrophoresis, proteins were transferred to a polyvinylidene fluoride (PVDF) membrane (Vazyme, China). The membrane was blocked, incubated overnight with a primary antibody against LILRA5 and β-actin (Proteintech, USA), and then incubated with a secondary antibody for 1 h. Results were detected using a Western Blot Imaging System (4000R, CareTream, USA).

Cecal ligation and puncture (CLP) models: Eight-week C57 mice were put into slumber after inhaling isoflurane, and a 22-gauge needle was used once to puncture a stump in order to release stool. The abdomen was then closed after the cecum was moved intraabdominally. Saline (0.2 mL) was injected via the abdomen for fluid resuscitation. Sham-operated mice were not ligated or punctured. Four hours following the procedure, the mice were put back in their cages and slaughtered. We induced sepsis in mice through intraperitoneal injection of LPS at a dose of 20 mg/kg, with some mice receiving PBS as a control group. This method established an LPS-induced sepsis model.

Intracellular ROS levels were quantified in THP-1 cells following LPS stimulation (10 ng/mL) or vehicle treatment (0.1% BSA). After a 24-h incubation at 37°C in a humidified 5% CO₂ incubator, ROS induction was assessed using the ROS assay kit (Beyotime, China), which employs the cell-permeant fluorogenic probe 2′-7′-dichlorofluorescein diacetate (DCFH-DA). This probe is hydrolyzed by cellular esterases into DCFH carboxylate, which is subsequently oxidized by intracellular ROS into the fluorescent 2′-7′-dichlorofluorescein (DCF). For the assay, cells were loaded with 15 μM DCFH-DA in culture media, incubated at 37°C for 30 min, and analyzed by flow cytometry without washing. The intracellular ROS levels were expressed as mean fluorescence intensity (MFI) values.

Data handling and visualization were performed using R 4.2.0. Statistical significance was determined using a two-tailed test, with p-values less than 0.05 considered significant.

Our study comprised a total of 21 samples, each demonstrating a uniform cell distribution. Given the consistency observed across all samples, we infer that batch effects likely had minimal influence on the results ( Supplementary Figure 1 ). Distributions of each sample on the single-cell level were visualized by the UMAP algorithm ( Figure 1A ). Based on single-cell analysis, cells were categorized into 23 different clusters ( Figure 1B ). Diverse cell types (monocytes, CD4⁺ T cells, B cells, neutrophils, NK cells, megakaryocytes, macrophages, CD8⁺ T cells, DCs, and mast cells) were unveiled according to different markers ( Figure 1C ). The densities of markers were shown (monocyte: S100A8, S100A12, CD14, and LYE; CD4⁺ T cell: CD3E, IL7R, CD27, and CCR7; B cell: CD79A and MS4A1; neutrophil: JAML and GOS2; NK cell: NKG7, GNLY, KLRB1, and KLRD1; megakaryocyte: PF4, GP9, and PPBP; macrophage: C1QA, C1QB, and CD68; CD8⁺ T cell: CD8A and CD8B; DC: FCER1A and CD1C; mast cell: GATA2, KIT, and CPA3) ( Figure 1D ). The density of each marker is shown in Supplementary Figure 2 . The features of markers are shown in Supplementary Figure 3 . Corresponding proportions of 10 cell types in different states of sepsis (health control, the early stage, and the late stage) were presented ( Figures 1E, F ). Gene ontology (GO) analysis revealed potential functional mechanisms and offered insights into the biological roles of 10 heterogeneous cell populations ( Supplementary Figure 4 ).

A total of 807 OS-related genes were downloaded from GeneCards with a relevance score ≥ 7. The OS was elevated in sepsis in both training and test sets ( Figures 2A, B ). Notably, the expression levels of the OS genes were significantly higher in both sepsis samples than in the control tissues. Supplementary Figure 5A shows the expression of 807 OS-related genes in three different stages of sepsis (healthy control, early stage, and late stage). The expression of the OS-related genes was significantly higher in the early stage of sepsis. Five algorithms (including AUCell, UCell, singscore, ssgsea, and AddModuleScore) were applied to evaluate the index of OS ( Figures 2C, D ). The expression of 807 OS-related genes in 10 types of cells is shown ( Figure 2E ; Supplementary Figure 5B ). According to the OS score, neutrophils, macrophages, DCs, and mast cells are four cell types with the most OS. Scores and density of OS-related genes in 10 types of cells are also visualized in Figures 2F, G . Moreover, OS scores of 10 cell types were also visualized in different stages of sepsis ( Figure 2H ). In contrast, neutrophils, macrophage, and DCs are three cell types with the most OS activity in the early stage of sepsis. In conclusion, neutrophils, macrophages, and DCs may play a potential role in response to sepsis.

To explore the cell type with the most OS, this study defined three cell groups by the quartile method according to the OS scores ( Figure 3A ). The top 25% cells were defined as the high-OS cells (25,963 cells with OS score > 2.419). The bottom 25% cells were defined as the low-OS cells (25,963 cells with OS score < 1.552). The medium cells were defined as those belonging to the “transition state”. Macrophages, neutrophils, DCs, and mast cells were four types of cells with the most OS in response to infection ( Figures 3B,C ). Moreover, we also define a new cluster of upregulation genes between the high-OS cells and the low-OS cells ( Figure 3D ). This new cluster of genes may have a real effect in those cells in response to sepsis.

We next applied hdWGCNA to further investigate the characteristics and functions of genes. After weighing median connectivity, mean connectivity, and scale-free topology model fit, a power value of 12 was selected ( Figure 3E ). Ten modules were generated accordingly ( Figure 3F ). Ten gene modules were obtained and the top huh gene was presented following the hdWGCNA pipeline ( Supplementary Figure 6A ). Correlations between every two modules are shown in Supplementary Figure 6B . Among 10 modules, the brown module showed the most positive correlation with macrophages and monocytes ( Figure 3G ). Genes in brown, yellow, turquoise, and purple modules showed a higher expression in high-OS cells ( Figure 3H ). Protein–protein interaction (PPI) networks of brown, yellow, turquoise, and purple modules were visualized ( Figure 3I ), while PPI networks of six other modules were also visualized ( Supplementary Figure 6C ). A total of 343 shared genes between OS-related genes and screened genes from hdWGCNA may be the real OS-related genes in sepsis ( Figure 3J ). Disease oncology (DO) and GO analysis of these hub genes were visualized ( Supplementary Figures 7A, B ).

To identify hub genes in the 343 shared genes, the differential expression analysis showed 292 DEGs ( Figure 4A ). The selected genes were then integrated into various machine learning algorithms, such as LASSO regression, SVM-REF, Boruta, GBM, and RF. We performed LASSO regression on these genes, which reduced the gene number to 12 ( Figure 4B ). The GBM algorithm screened 49 feature genes ( Figure 4C ). SVM-REF screened 26 feature genes. The accuracy and error of the SVM-FRE showed a good performance ( Figures 4D, E ). The RF algorithm screened 32 feature genes ( Figures 4F, G ). Moreover, the Boruta algorithms also found 77 hub genes ( Figures 4H, I ). Above all, the Venn diagram screened four shared genes (LILRA5, MGST1, PLBD1, and S100A9) as hub genes ( Figures 4J ).

To evaluate the reliability and accuracy, we evaluated the four feature genes (LILRA5, MGST1, PLBD1, and S100A9) in the bulk level. In the training dataset, the expression of LILRA5, MGST1, PLBD1, and S100A9 was upregulated significantly in the sepsis group compared to the normal group ( Figure 5A ). Moreover, receiver operating characteristic (ROC) curves were used to evaluate the diagnostic potential of the four genes. Area under the ROC curves for LILRA5, MGST1, PLBD1, and S100A9 were 0.980, 0.992, 0.995, and 0.998, respectively ( Figure 5B ). Moreover, in the test datasets, the expression of LILRA5, MGST1, PLBD1, and S100A9 was significantly elevated in the sepsis group ( Figure 5C ). In the test dataset, the areas under the ROC curve (AUCs) for LILRA5, MGST1, PLBD1, and S100A9 were 1.000, 0.979, 0.976, and 1.000, respectively ( Figure 5D ). To evaluate the correlation between four genes and OS, we utilized the bulk data for further analysis. Correlations between each of the two genes and OS-related genes in the training dataset and the test dataset were revealed ( Supplementary Figures 8A, B ). In the training and test dataset, OS activity was elevated in sepsis samples ( Supplementary Figures 8C, D ). These four genes were also significantly upregulated in the group with high OS activities ( Figure 5E ). Correlation of each of the four genes and the OS gene set in the training dataset and the test dataset revealed that these four genes were all positively associated with OS ( Supplementary Figures 8E, F ). Expression features of four genes in different cell types were also delineated at the single-cell level. The results demonstrated that S100A9, PLBD1, MGST1, and LILRA5 were highly expressed in DCs, macrophages, neutrophils, and monocytes ( Figure 5F ). Moreover, expressions of four genes in 10 types of cells were visualized ( Figure 5G ) and the results of UMAP also revealed the expression density of four genes ( Figure 5H ). Moreover, four genes were all positively associated with the OS activity at the single-cell level ( Supplementary Figure 7G ). In our study, we chose LILRA5 for further study. As LILRA5 was significantly upregulated in the macrophage, LILRA5 may play an important role in the OS activities in macrophages.

To explore the function of LILRA5 in macrophages, cell–cell communication and transcription factor analysis were performed. The overall communication ability of LILRA5⁺ macrophage and LILRA5⁻ macrophage with other cell types is shown in Figure 6A . Interactions between each cell type and other cell types are visualized in Supplementary Figure 9 . Interestingly, the outgoing interaction strength showed that the secretory ability of LILRA5⁺ macrophage was stronger compared with other cells (Figure 6B), including MIF, GALECTIN, and RESISTIN signaling (Figure 6C). Ligand–receptor analysis indicated that MIF-(CD74⁺CXCR4) and MIF-(CD74⁺CD44) were more activated in the signaling from LILRA5⁺ macrophage to B cells and DCs ( Supplementary Figure 10A ). Moreover, RETN-CAP1 signaling was most significant between neutrophils to LILRA5⁻ macrophage and neutrophils to LILRA5⁺ macrophage ( Supplementary Figure 10B ). Compared with LILRA5⁻ macrophage, LILRA5⁺ macrophage has a much higher proportion in the early stage of macrophage ( Figure 6D ). The pseudotime and the slingshot analysis showed that LILRA5⁺ was highly expressed in the early stage of macrophage ( Figures 6E, F ). Moreover, the CytoTRACE analysis was also performed. LILRA5⁺ macrophage had higher CytoTRACE scores, which indicated that this cell type appears in the earlier stage ( Figures 6G, H ). LILRA5⁺ macrophages have a higher tendency towards an “undifferentiated” state based on the predicted scores. Moreover, expression of LILRA5 correlated positively with cytoTRACE scores ( Figure 6I ). The results supported the conclusion that LILRA5 expression is higher in the early stages of development.

To validate the expression of LILRA5 in sepsis, we established the CLP mouse model, and the blood of mice was derived. Real-time quantitative reverse transcription (qRT-PCR) showed that LILIR5 was significantly upregulated in sepsis ( Figure 7A ). Our study proved that LILRA5 was mostly expressed in macrophages and potentially regulated OS activity. qRT-PCR showed that the mRNA expression of LILRA5 was higher following the stimulation of LPS in THP-1 cells ( Figure 7B ). The protein level of LILRA5 was higher ( Figure 7C ). To explore the role of LILRA5 in regulating OS in macrophage, we first silenced LILRA5 expression by siRNA ( Figure 7D ). After silencing the LILRA5 in THP-1 cells, the ROS levels were then measured. The result showed that the ROS level was lower when the LILRA5 was silenced ( Figure 7E ). In conclusion, LILIRA5 regulated OS of macrophages positively in sepsis.