The sepsis patient dataset GSE85233, comprising 22 normal and 51 sepsis samples, was retrieved from the Gene Expression Omnibus (GEO) database. Additionally, single-cell RNA sequencing data were obtained from the GEO dataset GSE167363. For independent validation, another dataset, GSE57065, including 25 control and 28 sepsis samples, was utilized.

A curated list of 203 coagulation-related genes was compiled using gene sets from the Gene Set Enrichment Analysis (GSEA) database. These gene sets were derived from the pathways hsa04610 and hsa04611, which are associated with coagulation and related processes.

The differential expression analysis of coagulation-related genes in sepsis patients was performed using the ‘limma’ package in R. To visualize the results, a volcano plot was generated with ‘ggpubr,’ and a heatmap depicting the expression levels of differentially expressed genes was created using the ‘heatmap’ package.

A correlation analysis examined the interrelationship among the 20 differently expressed genes. A correlation heatmap was generated in R using the “corrplot” package, and a circos plot created using the “RCircos” package was used to visualize the results (29).

Gene Ontology (GO) and Disease Ontology (DO) enrichment analysis was performed using the ‘clusterProfiler’ package to identify the biological functions and disease associations of the differentially expressed genes (30). Additionally, a Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was carried out using “clusterProfiler” to investigate the roles of genes in different biological pathways (31).

Using the “RCircos” package, which shows genomic data in a circular form and enables the identification of potential chromosomal patterns related to changes in gene expression in sepsis, the chromosomal positions of the 20 differentially expressed coagulation-related genes were determined.

A diagnostic model was constructed using machine learning techniques. The ‘randomForest’ package was used for random forest analysis to determine the key genes differentiating sepsis patients from controls. Support Vector Machine (SVM) analysis was conducted with the ‘SVM’ function within the ‘e1071’ package to optimize the accuracy and minimize error based on feature selection. The ‘glmnet’ package was utilized for lasso analysis to find diagnostic biomarkers by applying a penalty to the coefficient sizes. The common key genes among the different machine learning approaches were identified through a Venn diagram using the ‘Venn’ package.

The ‘limma’ package, intended for gene expression data analysis from microarray or RNA-seq technologies, validated the key gene expression. The ‘timeROC’ package was utilized to perform Receiver Operating Characteristic (ROC) curve analysis to evaluate the identified biomarkers’ diagnostic efficacy.

The ‘patchwork’ and ‘org.Hs.eg.db’ packages were used to conduct Gene Set Enrichment Analysis (GSEA) to investigate the key genes’ biological function and pathway enrichment analysis. Single-sample GSEA (ssGSEA) was used to identify important signaling pathways differently activated in sepsis patients compared to healthy controls. The analysis was done using ‘Limma’, and the correlation between the key genes and the signaling pathways was shown using ‘ggplot2’.

Following quality control filtering of single-cell RNA sequencing data from the GSE167363 dataset, the remaining cells were processed for dimensionality reduction and clustering. The ‘SingleR’ package was utilized for automated cell-type annotation, while ‘ggplot2’ was employed to assess the expression levels of key genes in different cell types.

The ‘CIBERSORT’ package, based on gene expression profiles and a predefined signature matrix of immune cell types, was used to deconvolute the immune cell proportions in sepsis patients. We visualized the differences in immune cell fractions between the control and sepsis samples using the ‘ggpubr’ package.

Based on the patterns of gene expression, different sepsis subtypes were distinguished using Principal Component Analysis (PCA) utilizing the ‘prcomp’ and ‘ggplot2’ functions. To classify sepsis patients into distinct subtypes, consensus clustering was carried out with the help of the ‘ConsensusClusterPlus’ package (32, 33). The ‘pheatmap’ package was used to create heatmaps of gene expression, and ‘ggpubr’ was used to generate box plots for visual comparison of differential gene expression across subtypes.

This study was approved by the Ethics Committee. Five sepsis patients and five healthy individuals undergoing routine health examinations were recruited at Nanfang Hospital, Southern Medical University, between December 2023 and January 2024. Peripheral blood mononuclear cells (PBMCs) were isolated from patients’ peripheral blood samples using previously described methods (34). Total RNA was extracted from PBMC samples using the FastPure Cell/Tissue Total RNA Isolation Kit (Vazyme). RNA was then reverse-transcribed into cDNA using the ReverTra Ace qPCR RT Master Mix and gDNA Remover Kit. Quantitative real-time PCR (qRT-PCR) was performed using the SYBR Premix Ex Taq II in a real-time fluorescence quantitative PCR system, with GAPDH selected as the endogenous control for mRNA. The reaction conditions were as follows: initial denaturation at 95°C for 10 minutes, followed by 45 cycles of 95°C for 5 seconds and 60°C for 30 seconds (35). The amplification of target genes and internal reference genes was performed separately for each sample, with each group of samples containing three replicate wells. Data analysis was conducted using the 2^(-ΔΔCt) method. The primer sequences are provided in Supplementary File 1 .

Rstudio was used for all statistical analysis and computational modeling (version 4.2.1).

We studied 203 coagulation genes using differential gene expression analysis to comprehend the connection between coagulation-related genes and sepsis. This analysis revealed 20 genes that exhibited significant differential expression, as shown in Figures 1A, B . Furthermore, we conducted a correlation study to explore the relationships among these genes. The results from the circos plot indicated a substantial correlation among all 20 genes, as depicted in Figure 1C . Additionally, the heatmap of correlations reinforced these findings, demonstrating a close association between these 20 genes, as illustrated in Figure 1D .

To elucidate the biological roles of the differentially expressed genes, we conducted Gene Ontology (GO) analysis and pathway enrichment analysis using the Kyoto Encyclopedia of Genes and Genomes (KEGG). The GO analysis indicated that these genes are primarily involved in biological functions such as immune receptor activity, coagulation, secretory granule lumen composition, and immune response activation ( Figure 2A ). The Disease Ontology (DO) enrichment analysis revealed that these genes are enriched in pathways related to tumor signaling and autoimmune hemolytic anemia ( Figure 2B ). KEGG pathway analysis demonstrated that these differentially expressed genes are mainly concentrated in pathways associated with complement system activation and platelet activation ( Figure 2C ). Additionally, the chromosomal locations of these 20 genes were visualized using a circos plot ( Figure 2D ).

Our results indicated a significant relationship between sepsis and immune responses. Further investigation through CIBERSORT analysis elucidated the differences in immune cell composition in sepsis patients. The findings revealed that sepsis patients exhibited lower levels of unactivated B cells, CD8 T cells, and unactivated NK killer cells while having higher proportions of monocytes and neutrophils, as shown in Figures 3A, B . Subsequent correlation analysis between the 20 differentially expressed genes and immune cells demonstrated that most genes significantly associate with immune cells, as depicted in Figure 3C .

We employed consensus clustering through differential gene expression analysis to categorize sepsis patients into two subtypes, as illustrated in Figures 4A–C . Analysis of the differential genes between these two subtypes revealed that most of the 20 genes exhibited significant differences, as presented in Figures 4D, E .

Principal Component Analysis (PCA) results suggested a good distinction between the two sepsis subtypes, as indicated in Figure 5A . An examination of immune cell infiltration in these subtypes revealed that subtype A had higher immune cell infiltration, as shown in Figure 5B . A heatmap of correlations highlighted significant relationships between immune cells and the 20 genes, as seen in Figure 5C .

To create a diagnostic model, we used machine learning techniques. As shown in Figures 6A, B , random forest analysis identified the top five key genes, including FCER1G and A2M. The Support Vector Machine (SVM) results indicated that the selection of two genes offered the best accuracy and minimal error. However, we chose a set of five genes for further analysis, achieving an accuracy of 1 and an error rate of 0, as depicted in Figures 6C, D . Lasso analysis identified six diagnostic biomarkers, as shown in Figure 6E . To develop a common key gene diagnostic model, a Venn diagram was utilized to find the intersecting genes among the three analyses, resulting in two intersecting genes: FYN and FCER1G, as illustrated in Figure 6F . The qRT-PCR results showed that FYN was downregulated in the sepsis patient group compared to the healthy control group ( Figure 7A ), while FCER1G was upregulated in the sepsis patient group compared to the healthy control group ( Figure 7B ).

We assessed the diagnostic capability of the two identified biomarkers through column line graphs, as seen in Figure 8A . Decision Curve Analysis (DCA) suggested that patients could benefit from these biomarkers, as shown in Figure 8B . Calibration curves indicated a minimal difference between the actual risk of sepsis and the predicted risk, signifying the high accuracy of the model, as presented in Figure 8C . The ROC curves of the column line graphs demonstrated good predictive performance, as depicted in Figure 8D . Additionally, the diagnostic ROC for the key genes, FYN and FCER1G, showed areas under the curve (AUC) of 0.999 and 1.000, respectively, indicating excellent diagnostic efficacy, as shown in Figures 8E, F .

GSEA enrichment analysis revealed that FCER1G was predominantly associated with pathways such as fatty acid synthesis and the intestinal immune network for IgA production ( Figure 9A ). In contrast, FYN was primarily enriched in pathways related to transplant rejection, antigen processing, and presentation ( Figure 9B ). Using ssGSEA, we examined differences between sepsis patients and healthy controls across 50 hallmark signaling pathways. In sepsis patients, several hallmark pathways, including KRAS_SIGNALING_DNE, REACTIVE_OXYGEN_SPECIES_PATHWAY, and TNFA_SIGNALING_VIA_NFKB, were significantly upregulated ( Figure 9C ). Additionally, a correlation analysis between the two key genes and various hallmark pathways revealed significant associations with most pathways ( Figure 9D ).

The expression levels of the key genes and their diagnostic efficacy were validated in a separate cohort. FCER1G expression was found to be elevated in sepsis patients ( Figure 10A ), while FYN expression was reduced ( Figure 10B ), consistent with findings from the training set. ROC curve analysis demonstrated the strong diagnostic potential of these genes, with the area under the curve (AUC) being 1 for FCER1G and 0.985 for FYN ( Figures 10C, D ).

After quality control filtering, 25,458 cells were extracted from the GSE167363 dataset. Post-dimensionality reduction and clustering, these cells were categorized into 11 distinct clusters ( Figure 11A ). Further cell annotation identified six cell types ( Figure 11B ). FCER1G was found to be highly expressed in monocytes, NK cells, and platelets ( Figure 11C ), whereas FYN showed higher expression in CD4+ T cells and NK cells ( Figure 11D ).