We obtained 181 glycosylation-related genes from the GlycoGene DataBase (GGDB) database (https://acgg.asia/ggdb2/). We downloaded sepsis-related datasets from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) (16): the training cohort GSE65682, which consists of transcriptome sequencing data and corresponding follow-up data from whole blood samples isolated from 42 healthy controls and 760 sepsis subjects. Additionally, the datasets GSE54514 and GSE95233, which include follow-up information, were used to validate the prognostic value of the predictive model. Furthermore, six other sepsis-related datasets (GSE131761, GSE137340, GSE236713, GSE28750, GSE570653, GSE69528) were chose to confirm the diagnostic value of the model. To further investigate the underlying mechanisms, we downloaded GSE1754538 from GEO, which contains single-cell RNA sequencing (scRNA-seq) data from whole blood samples collected from 5 healthy controls and 4 sepsis donors. To validate our findings, GSE176363 was enrolled, which contains 2 healthy controls and 10 sepsis donors. The scRNA-seq data was processed using Seurat for data loading, quality control, and dimensionality reduction. Cell type annotation was performed using the “SingleR” package and cell specific markers were obtained via the FindAllMarkers function to validate the cell type annotation results.

We prospectively enrolled 30 patients diagnosed with sepsis and 20 healthy controls as validation dataset. Patients were recruited from the intensive care units of Peking University People’s Hospital from January 2023 to April 2024. The inclusion criteria were based on the Sepsis-3.0 definitions, which require a suspected infection plus a Sequential Organ Failure Assessment (SOFA) score of 2 or higher. Healthy controls were selected from individuals undergoing routine health check-ups who did not have any signs or symptoms of infection or inflammation.

We used the “limma” package to identify genes differentially expressed between normal and sepsis patient blood samples (P < 0.05 & |LogFC| > 1) (17). From these, we selected glycosylation-related genes for further analysis. These selected genes underwent functional enrichment analysis using the Metascape database (www.metascape.org/) (18), with P < 0.05 considered significant. Finally, we used univariable Cox regression analysis to identify glycosylation-related genes significantly affecting sepsis prognosis, designating them as Hub-Glys.

We used the “ConsensusClusterPlus” package in R to perform consensus clustering based on Hub-Glys (19). The optimal number of clusters was determined by analyzing the cumulative distribution function (CDF) curve. Subsequently, we used the “prcomp” function in R for principal component analysis (PCA) to confirm the reliability of the clustering. Kaplan-Meier (KM) survival analysis and log-rank tests were conducted using the survival package in R. The survival differences between different groups were analyzed, and the survival status of all patient samples within 25 days was displayed. Immediately after, we identified glycosylation-related differentially expressed genes (Gly-DEGs) through differential analysis between diverse sepsis groups. These candidate genes were then subjected to gene set enrichment analysis (GSEA) using the clusterProfiler 4.8.2 package in R (20), with the Reactome pathway database for functional enrichment. A P value < 0.05 was considered significantly enriched.

Based on GSE65682 cohort, we used three machine learning algorithms to select the prognostic Gly-DEGs. The “randomForestSRC” package in R was obtained for random survival forest (RSF) analysis (21). An ensemble of 1000 survival trees with default parameters was conducted and the out-of-bag (OOB) error was used to assess model stability. Variables were ranked by minimal depth and variable importance (VIMP). Genes with positive VIMP and consistent presence across trees were retained as candidate prognostic features. We performed CoxBoost algorithm analysis (22), via the CoxBoost package with 10-fold cross-validation to determine the optimal number of boosting steps. Genes with non-zero coefficients were selected. The “glmnet” package for The LASSO regression analysis (23) was conducted using 10-fold cross-validation to determine the optimal λ. Genes with non-zero coefficients at λmin were retained. To select genes for the final prognostic model, we applied LASSO, CoxBoost, and RSF to the Gly-DEGs. Selected genes filtered by at least two of the three algorithms were considered robust predictors and included in the final multivariable Cox regression analysis to establish the Sepsis Risk Score (GRS). This integrative strategy reduced model overfitting and emphasized features supported across multiple selection methods. Based on the GRS, patients were stratified into high and low-risk groups according to the median GRS value, and their prognostic differences were analyzed in GSE65682 and validated in the GSE54514 as well as GSE95233 datasets.

To enhance the clarity of the diagnostic and prognostic value of the GRS, we utilized the “meta” package in R to combine the Odds Ratios (OR) and Hazard Ratios (HR) obtained from multiple studies. Depending on the degree of heterogeneity determined by the I² statistics, we extracted and pooled the data from each study using either a fixed-effect or a random-effects model.

We utilized the heatmap package to visualize the distribution of clinical characteristics within the GRS groups from the GSE65682 and GSE95233 datasets. Additionally, we employed univariable and multivariable Cox regression to select independent prognostic factors, thereby assessing the independent prognostic performance of the GRS. Furthermore, we evaluated the diagnostic capability of the GRS and its model genes for sepsis by using the “pROC” package to plot Receiver Operating Characteristic (ROC) curves across eight datasets (24).

We used the “Seurat” R package to preprocess and analyze scRNA-seq data. Normalization was performed using the “NormalizeData” function with the “LogNormalize” method, and the data was converted into a Seurat object. For quality control, we calculated the percentages of mitochondrial and ribosomal genes, excluding cells with fewer than 200 or more than 3000 genes, and those with over 20% ribosomal RNA content. “Harmony” function was applied to correct for potential batch effects. Since all scRNA-seq samples were derived from the same GEO dataset (GSE175453) with unified experimental conditions, no significant batch effect was observed based on UMAP distribution, and batch correction was not applied. We selected the top 3000 variable features using the “FindVariableFeatures” function and standardized the scRNA-seq data. The “ScaleData” and “RunPCA” functions were used to obtain principal components (PCs), followed by dimensionality reduction using Uniform Manifold Approximation and Projection (UMAP). Cell types and subtypes were automatically annotated using the “SingleR” package with the HumanPrimaryCellAtlasData as reference, then further validated using cell-type–specific differential gene expression analysis with FindAllMarkers function. Based on these DEGs and their functional relevance, we confidently validated the major cell types. This data-driven approach provides robust support for cell identity assignment and avoids potential circularity associated with relying solely on predefined marker lists.

Total RNA was extracted using the TRIzol™ Plus RNA purification kit and reverse transcribed using a reverse transcription kit. RT-QPCR was performed on a AriaMx Real-time PCR System using SYBR^® Green real-time PCR master mix. Gene expression was analyzed using the 2ΔΔCt method with the threshold cycle (Ct). Data were normalized to the expression of GAPDH, and reported as expression in controls set to 1. The primer sequences are listed in the Supplementary Table 1 .

All statistical analyses were conducted using R software environment (version 4.3.1). The Wilcoxon rank-sum test was used to compare differences between two groups. The Kruskal-Wallis test was employed to evaluate differences among more than two groups. Spearman’s correlation was utilized for correlation analysis. A P-value of less than 0.05 was considered statistically significant. The significance levels were denoted as *p < 0.05; **p < 0.01; ***p < 0.001.

The workflow of the present study is illustrated in Figure 1 . To identify potential glycosylation-related biomarkers in sepsis, we first screened for differentially expressed genes using the “limma” R package on the GSE65682 dataset, which contains RNA-seq data from peripheral blood of sepsis patients and normal controls. This identified 2,774 differentially expressed genes, with 952 being upregulated and 1,822 being downregulated in sepsis patients ( Figure 2A ). We further screened for glycosylation-related genes based on the GGDB database, identifying 38 differentially expressed genes ( Figure 2B ). Among these, 23 were upregulated and 15 were downregulated in sepsis patients ( Figure 2C ). Functional enrichment analysis using the Metascape database revealed that these 38 genes were significantly enriched in multiple glycosylation-related pathways ( Figure 2D ).

To identify Hub-Gly affecting patient prognosis, we conducted a univariable Cox regression analysis, selecting 8 Hub-Glys with significant prognostic impact ( Figure 3A ). Based on these Hub-Glys, consensus clustering analysis was performed, dividing the 802 samples in GSE65682 into two subtypes (Subtype1 and Subtype2) ( Figure 3B ). PCA showed clear discrimination between these subtypes ( Figure 3C ). Disease status distribution indicated that healthy patients were concentrated in Subtype1, with a significant difference in sepsis status between the subtypes (p=1.15e-30) ( Figure 3D ). After excluding healthy and duplicate samples, Kaplan-Meier survival analysis was conducted on 479 sepsis patients. The survival curves demonstrated a significant difference in prognosis, with Subtype1 patients exhibiting better outcomes ( Figure 3E ).

We conducted differential expression analysis between the two subtypes and identified 100 Gly-DEGs ( Figure 4A ). Reactome pathway enrichment analysis revealed that these Gly-DEGs are significantly enriched in glycosylation-related pathways, such as REACTOME O-Linked Glycosylation ( Figure 4B ) and REACTOME O-Linked Glycosylation of Mucins ( Figure 4C ).

Using the training set GSE65682, we applied multiple machine learning algorithms to select key glycosylation-related prognostic factors from these 100 Gly-DEGs. The Lasso algorithm identified 18 important prognostic factors (The full names of all gene abbreviations in this article can be found in Supplementary Table 2 ): MAN1A1, EXT2, GALNT11, CD24, EXT1, MOGS, HMMR, HIF1A, CD44, SELL, MUC1, SDC4, PDGFB, SIGLEC7, GALNT10, GCNT2, HBEGF, TNF ( Figure 4D ). The CoxBoost algorithm identified 20 important factors: MAN1A1, RPS6, EXT2, GALNT11, GALNT14, EXT1, MOGS, HMMR, FUCA1, AKT2, HIF1A, CD44, SELL, MUC1, SDC4, PDGFB, SIGLEC7, GCNT2, HBEGF, TNF ( Figure 4E ). The RSF algorithm identified 4 important factors: MAN1A1, SELL, HMMR, TNFRSF25 ( Figure 4F ). Integrating results from multiple algorithms enhances the stability of the prognostic model. Thus, we considered glycosylation-related prognostic factors selected by at least two algorithms as potential clinically relevant prognostic factors, yielding a total of 16 potential prognostic factors for sepsis ( Figure 4G ). Multivariable Cox regression analysis of these factors established the following sepsis prognosis model based on their coefficients:

GRS = CD44*(-0.46) + EXT1*10.52 + EXT2* (-0.26) + GALNT11*(-0.16) + HIF1A*0.50 + HMMR*0.23 + MUC*10.16 + SDC*40.20 + SELL*(-0.42),

where 6 factors demonstrated independent prognostic ability: CD44, EXT1, EXT2, HIF1A, HMMR, SELL ( Figure 4H ). Their functions in immune regulation or sepsis reported by previous literatures were summarized in Supplementary Table 3 .

To verify the reliability of the GRS in predicting patient prognosis, we conducted external validation using the datasets of GSE65682, GSE95233 and GSE54514. Sepsis patients were divided into high GRS and low GRS risk groups based on the median GRS value in the training set GSE65682 ( Figure 5A ). Significant survival differences were observed between the high and low GRS groups in GSE65682 (p<0.0001) ( Figure 5B ), GSE95233 (P=0.011) ( Figure 5C ), and GSE54514 (p=0.0098) ( Figure 5D ). The bar chart revealed significant distribution differences in survival states between the high and low GRS groups. In the GSE65682 dataset, High GRS group patients accounted for 75% of fatal cases, with a significant difference between groups (p=1.27e-08) ( Figure 5E ). In the GSE95233 dataset, High GRS group patients accounted for 76% of fatal cases, with a significant difference (p=0.04) ( Figure 5F ). Finally, in the GSE54514 dataset, High GRS group patients accounted for 80% of fatal cases, with a significant difference (p=0.01) ( Figure 5G ).

To explore the distribution of clinical characteristics among patients with different GRS states, we created a heatmap of clinical information for dataset GSE65682 ( Supplementary Figure S1A ) and dataset GSE95233 ( Supplementary Figure S1B ). The results showed that GRS was significantly associated with ICU-acquired infection (p<0.01) and gender (p<0.05). Furthermore, to determine the independent prognostic ability of GRS, we performed univariable Cox regression analysis ( Supplementary Figure S1C ) and multivariable Cox regression analysis ( Supplementary Figure S1D ) in GSE65682, incorporating GRS and clinical characteristics. The results indicated that GRS could serve as an independent prognostic factor for sepsis patients. This finding was further validated externally in GSE95233 ( Supplementary Figures S1E, F ), confirming its independent prognostic capability. To evaluate the generalizability and discriminative power of predictive factors across datasets, we calculated the C-index for each model in these three independent sepsis-related cohorts. As shown in Supplementary Figure S3A , C-index values varied across both datasets and predictors, reflecting differences in model robustness and cohort characteristics. GRS achieved high C-index values across multiple datasets, indicating strong and stable discriminatory performance.

To further investigate the role of GRS in diagnosing sepsis, we generated ROC curves using GRS across multiple datasets containing both healthy individuals and sepsis patients. The results showed that the AUC values for GRS were consistently around 0.80, indicating excellent diagnostic capability ( Figures 6A–H ). This suggests that GRS is a reliable tool for diagnosing sepsis. Additionally, a meta-analysis revealed that the summary OR for GRS was 24.89 (95% CI = 2.66-232.95) ( Figure 6I ). We also assessed the diagnostic performance of six independent prognostic factors. Supplementary Figures S2A-F show that CD44, EXT1, EXT2, HIF1A, HMMR, and SELL all demonstrated good diagnostic performance.

To reveal glycosylation differences between sepsis patients and healthy individuals at the single-cell level, we analyzed 21,644 cell samples from four sepsis patients and 27,808 cell samples from five healthy controls. No obvious batch effects were observed in tSNE plot of different samples ( Supplementary Figure S3B ). We identified seven cell types in both normal and sepsis samples: hematopoietic stem cells (HSCs), monocytes, platelets, T cells, natural killer (NK) cells, B cells, and neutrophils ( Figures 7A, B ). Cell-type–specific differential gene expression analysis validated the annotation results ( Supplementary Table 4 ).

We assessed the expression of six independent prognostic factors in these cell types and found varying levels of expression in all cell types except platelets ( Figures 7C, D ). Differential expression analysis showed that, except for CD44, the other five genes exhibited significant differences in monocytes ( Figure 7E ). All six genes showed significant differences in T cells ( Figure 7F ). In NK cells, HIF1A, HMMR, and SELL showed significant differences ( Figure 7G ). In B cells, CD44, EXT2, and HMMR were significantly different ( Figure 7H ). In HSCs, all genes except EXT2 and HMMR showed significant differences ( Figure 7I ). In neutrophils, HIF1A and SELL showed significant differences ( Figure 7J ). The details of Wilcoxon rank-sum test results were presented on Supplementary Table 5 . To validate our findings, we applied the same analysis in GSE167363 ( Supplementary Figures S3C, D , S4 ). Significant differences were also observed. Overall, the six independent prognostic factors showed significant differential expression in most major cell types, which may underlie their strong diagnostic capability. GSEA enrichment analysis revealed significant differences in the N-linked glycosylation pathway between normal and sepsis samples ( Figure 7K ).

Aiming at the glycosylation-related immune dysfunction, we visualized the differentially expressed glycosylation-related genes and performed enrichment analysis among immune cell types ( Figure 8A ). Monocytes exhibited high expression of TIMP1, VIM, LGALS3, and LGALS1, enriched in pathways related to immune activation and antigen presentation, indicating glycosylation’s role in monocyte-driven inflammatory responses. T cells upregulated MYC and RPS6, enriched in cell cycle and lysosomal pathways, suggesting glycosylation’s involvement in T cell proliferation and exhaustion priming. NK cells showed elevated expression of GNPTAB and KDELR2, with enrichment in N-glycan processing and pathogen response pathways, linking glycosylation to cytotoxic granule trafficking and effector function. B cells demonstrated ER-related enrichment in protein folding and retrograde transport, implicating glycosylation in antibody maturation. These findings reveal cell-type-specific glycosylation signatures that contribute to immune dysfunction in sepsis. Genes such as MYC (a transcription factor, TF), SELL, and LGALS3 (associated with cold shock domain-containing protein A) may act as upstream regulators or potential therapeutic targets in the context of glycosylation dysfunction related to sepsis. Then, we performed correlation analysis between six GRS genes and exhaustion markers in T cells ( Figures 8B, C ). We observed stronger correlations between glycosylation-related genes and T cell exhaustion markers in normal T cells compared to sepsis-derived T cells. This suggests that glycosylation may normally act in coordination with exhaustion-related transcriptional programs, but this relationship is disrupted under septic conditions, possibly contributing to dysfunctional immune responses.

We enrolled 30 sepsis patients, with their clinical information listed in Supplementary Table 6 . We validated the expression of key glycosylation-related genes in peripheral blood cells through PCR. The results showed that, compared with the healthy control group, the expression of EXT1, HIF1A and HMMR was upregulated in sepsis, whereas the expression of CD44, EXT2 and SELL was downregulated in sepsis group ( Figures 9A–F ).