We obtained a sepsis RNA microarray dataset and a sepsis RNA-seq dataset [GSE28750, GSE154918], respectively (Sutherland et al., 2011; Herwanto et al., 2021). Among them, there were 30 samples (20 normal and 10 sepsis) in the GSE28750 dataset. In the GSE154918 dataset, 60 samples (40 normal and 20 sepsis) were selected. Two samples of peripheral blood mononuclear cells (PBMCs) with sepsis and two samples of normal PBMCs were selected in the GSE167363 dataset (Qiu et al., 2021) from GEO database (Barrett et al., 2013); we also obtained the MetS dataset GSE98895 (D'Amore et al., 2018). Forty samples (20 normal and 20 MetS) were selected in the GSE98895 dataset. In order to analyze the data, the software R was employed. Principal component analysis (PCA) was used before and after correcting batch effects and visualizing the distribution of these datasets. The datasets were corrected for background, transformed by log2, and normalized. We also merged the datasets GSE28750 and GSE154918, and used the Combat method in the “sva” package to batch-correct the merged data (Leek et al., 2012). Then, the merged result was viewed by PCA dimensionality reduction algorithm. Following the identification of differentially expressed genes (DEGs) by the limma analysis, the following filtering criteria was applied to screen for significantly differentiating genes: p < 0.05 and |log2 Fold change (FC) | > 1.0 (Ritchie et al., 2015). Finally, Venn plots were used to illustrate common genes which are up- and downregulated in Sepsis and MetS, respectively. Meanwhile, Gene Ontology (GO) enrichment analysis of the common DEGs were conducted utilizing the “clusterProfiler” package of R software (Wu T. et al., 2021). The pathways with p < 0.05 were considered significant.

The “WGCNA” R package was employed to identify genes associated with Sepsis and MetS using a weighted gene co-expression network (Langfelder and Horvath, 2008). We used genes with expression >0 for further analysis to exclude outlier data. For the construction of the co-expression network, co-expression analysis was utilized. Flash clust was used for cluster analysis in our study. Clustering each sample from the beginning ensured the reliability of the network. As a result of calculating Pearson product moment correlation coefficients between gene pairs, we group genes with similar expression patterns into modules, thus creating a correlation matrix. Soft threshold functions are also used to transform the correlation matrix into a weighted adjacency matrix. To identify the most relevant Sepsis and MetS modules, we set the optimal soft threshold and identified the multi-co-expressed module genes simultaneously.

In order to further filter candidate genes for Sepsis and MetS diagnosis, two machine learning algorithms have been applied: random forest (RF) (Garge et al., 2013) and least absolute shrinkage selection (LASSO) (Alhamzawi and Ali, 2018). The search for important genes was carried out using the “random forest” R package. Based on decision tree theory, the RF algorithm was classified according to its ability to handle high-dimensional data and select highly informative gene clusters. The RF algorithm was used to screen diagnostic genes whose importance score exceeded 0.5. LASSO regression can be used for high-dimensional data to enhance the effectiveness of the analysis. For further reducing the dimensions of the obtained genes, the LASSO algorithm was applied to obtain the final diagnostic genes. The genes deemed most significant were selected as the core genes for further research.

In supervised integrated learning, eXtreme Gradient Boosting (XGBoost) is one of the most commonly used algorithms due to its scalability and convenience (Parente, 2021). For the model, hyperparameters were tuned using an optimisation method based on a Bayesian sequence model. Optimisation was performed on the training set, using K-fold cross-validation (K = 10) for continuous iteration. Candidate gene models were constructed using XGBoost on the training data set (GSE154918) and evaluated on the validation data set (GSE28750). Following that, the diagnostic efficacy of the model was evaluated using receiver operating characteristic curves (ROC), precision-recall curves (PR), and areas under the curve (AUC). This was verified in MetS patients.

Assess the presence of immune cells in each sample using the CIBERSORT algorithm (Hu et al., 2022). Based on linear support vector regression, the CIBERSORT deconvolution algorithm calculates the percentage of 22 immune cells in tissues or cells using machine learning. These 22 cell types included dendritic cells, CD4⁺ and CD8⁺ T Cells, B cells, macrophages M1 and M2, monocytes, neutrophils, natural killer cells, and natural killer cells. The proportion of immune cells in peripheral blood mononuclear cells (PBMCs) was compared between the disease and control groups. Meanwhile, the relationship between the Hub gene and immune cells in Sepsis and MetS was explored.

To analyze the dataset of single cells GSE167363, we used the “Seurat” R package to run PCA and t-distributed stochastic neighbor embedding (t-SNE) (Butler et al., 2018). Those cells with more than 4,000 features, mitochondrial genes over 25%, or less than 200 features were excluded from the analysis. After scaling the level of gene expression, a technique called “LogNormalize” was utilized to normalize the data. After normalizing the data, 3,000 highly variable genes (HVGs) within each sample were identified using the “vst” method. Principal component analysis (PCA) was then performed and the significant principal components (PCs) were identified using the elbow method. In the end, t-SNE analysis was performed using 20 PCs that were chosen. We used FindClusters function to cluster cells into 21 clusters. In order to locate differentially expressed genes (DEGs) for each cluster, the logfc. threshold parameter was set to 0.25 using the FindeMarker function. The “singleR” package comes with seven reference datasets, of which 5 are human and 2 are mouse, and we have selected “HumanPrimaryCellAtlasData” as the reference dataset (Aran et al., 2019). An automated annotation using “SingleR”package in conjunction with DEGs in each cluster to identify cell types, and then identifying the cell types in each cluster separately (Li et al., 2020). The Hub gene expression was also visualized by violin diagrams in different immune cells.

Molecular Signature Database (MSigDB) (Li J. et al., 2021) was used to download the hallmark gene set, and single sample gene set enrichment analysis (ssGSEA) was done to analyze metabolic pathways associated with the Hub genes. Lastly, we analyzed the correlation between immune cells and metabolism in Sepsis and MetS. Use the Pearson correlation method in the “stats” package of R language to calculate its correlation. Differences in metabolic pathway scores between single cell subpopulations were demonstrated by violin plots, where significant differences were determined by Wilcoxon tests.

All statistical tests were performed using R version 4.1.2. The Wilcoxon or Student’s t-test was used to analyse the difference between the two groups. Correlations between variables were determined using Pearson’s or Spearman’s correlation test. Statistical significance was set at a two-tailed p < 0.05.

As shown in Figure 1, the study flow chart explains how it was conducted. The PCA was performed on three datasets (GSE28750, GSE154918 and GSE98895) before corrections and normalizations (Supplementary Figures S1A, B). The datasets were normalized, and 3,902 DEGs (1930 upregulated and 1972 downregulated) were found in Sepsis, while 2,639 DEGs (1,354 upregulated and 1,285 downregulated) were found in MetS. By identifying common DEGs between Sepsis and MetS, 122 common upregulated DEGs and 90 common downregulated DEGs were found (Figures 2A, B) (Supplementary Tables S1, S2). GO enrichment analysis of the identified common DEGs was performed to investigate their biological functions and pathways. According to GO analysis, the commonly upregulated DEGs are mainly involved in cell activation and leukocyte activation involved in immune response and regulation of regulated secretory pathway, while the common downregulated DEGs were enriched in epithelial tube morphogenesis, actin cytoskeleton, mitochondrial matrix, SMAD protein signal transduction (Figures 2C, D) (Supplementary Tables S3,S4).

With a threshold of 80, 2 outlier samples were detected and removed, and 98 samples were retained (Supplementary Figures S1C, D). The “pick Soft Threshold” function of the “WGCNA” package is used to filter out power parameters from 1 to 30. As a soft threshold, a power of 6 was selected for ensuring the scale-free network (Figure 3A). A total of 14 modules containing genes with similar co-expression traits were obtained using the “cutree” dynamic and module eigengenes functions (Figure 3B). The heatmap displayed the correlation between each module and the diseases (Figure 3C). “Brown” modules indicate that Sepsis and MetS are highly linked (Sepsis: r = 0.46, p = 0.009; MetS: r = 0.26, p = 0.003). Sepsis and MetS have positively linked genes in the brown module (Sepsis: cor = 0.38, p = 2.8e-18; MetS: cor = 0.37, p = 2.4e-17) (Figures 3D,E). For this brown module gene, a GO analysis was performed. The results show that it was mainly enriched in histone modification, peptidyl−lysine modification, regulation of response to DNA damage stimulus in biological process (BP), Mitochondrial matrix, mitochondrial inner membrane and Mitochondria containing protein complexes in cellular component (CC) and transcription coregulator activity and structural constituent of ribosome in molecular function (MF) (Figure 3F).

We used the RF algorithm in combination with LASSO regression to finally obtain seven diagnostic genes, including STOM, BATF, CASP4, MAP3K14, MT1F, CFLAR, UROD (Figures 4A–D). Afterwards, we evaluated these genes’ diagnostic value. The AUC values of ROC curves were 0.995 of STOM (Supplementary Figure S2A), 0.996 of BATF (Supplementary Figure S2B),0.995 of CASP4 (Supplementary Figure S2C), 0.995 of MAP3K14 (Supplementary Figure S2D), 0.968 of MT1F (Supplementary Figure S2E), 0.934 of CFLAR (Supplementary Figure S2E), 0.973 of UROD (Supplementary Figure S2E). All seven gene features had high accuracy with AUC >0.9, demonstrating their predictive power. Based on the training set GSE154918, we constructed a candidate gene model (STOM, BATF, CASP4, MAP3K14, MT1F, CFLAR) and evaluated it on the validation set GSE28750. As displayed in Figure 4E, in GSE154918, the AUC of ROC value was 0.997 and the PR value was 0.995. The ROC and PR values for GSE28750 are 0.965 and 0.951, respectively (Figure 4F), demonstrating the model’s diagnostic accuracy. It has also been validated in MetS, indicating that the model is applicable and effective in MetS, with a ROC of 0.97 and a PR of 0.971 (Figure 4G).

Sepsis and MetS patients with immune infiltration were studied. In addition, heat maps show the differential expression of seven key genes in immune cells (Figures 5B, D). Normal tissues contained fewer neutrophils and monocytes than Sepsis tissues (p < 0.05). A comparison of Sepsis tissues and normal tissues revealed that Sepsis tissues contained significantly fewer naïve B cells, memory naïve B cells, CD8 naïve T Cells, and CD4 naïve T Cells (Figure 5A). The expression of STOM, BATF, CASP4, MT1F, CFLAR, and UROD was negatively correlated with infiltration levels of resting NK cells, CD4 naïve T Cells, CD8 T Cells, and CD4 resting T Cells. And MAP3K14 expression was negatively associated with neutrophils, activated mast cells, monocytes, macrophage M0, and NK activated cells (Figure 5B). We also calculated the immune cell content in patients with MetS, and monocyte proportions were higher in patients with MetS than in controls (Figure 5C). The expression levels of key genes differed significantly in patients with MetS, with six genes (STOM, BATF, CASP4, MT1F, CFLAR, UROD) being expressed at lower levels in M0, M1 and M2 macrophages than in MAP3K14. However, six genes (STOM, BATF, CASP4, MT1F, CFLAR, UROD) were expressed at higher levels in monocytes, B memory cells and T cell regulation (Tregs), but all higher than MAP3K14 (p < 0.05) (Figure 5D). To explore the potential metabolic pathways involved in hub genes, correlations between hub genes and classical metabolic pathways were analysed. In sepsis samples, significant positive correlations were found between six hub genes (STOM, BATF, CASP4, MT1F, CFLAR, UROD) and the pathways of glycolysis, bile acid metabolism, adipogenesis, cholesterol homeostasis and xenobiotic metabolism, while the pathways of glycolysis, bile acid metabolism, adipogenesis, cholesterol homeostasis and xenobiotic metabolism showed negative correlations with MAP3K14 (Figure 5E). Similarly, MetS samples expressing six central genes (STOM, BATF, CASP4, MT1F, CFLAR, UROD) showed positive correlations in glycolysis, oxidative phosphorylation and fatty acid metabolism pathways (Figure 5F), while MAP3K14 expression was negatively correlated with glycolysis, oxidative phosphorylation and fatty acid metabolism (Figure 5F).

In order to check the quality of the single-cell dataset GSE167363, a preliminary quality check was performed. The correlation between nFeature RNA, nCount RNA, and precent. mt was examined to make sure the cell samples used in the study were of high quality. Figure 6A exhibited a positive correlation between nCount RNA and nFeature RNA representing unique molecular identifiers, with a correlation coefficient of 0.94. We excluded some cells and the result was diaplayed in Figures 6B,C. In the scRNA-seq dataset, a total of 3,000 genes with high levels of variation were identified, and ten of the most significant markers were tagged. A PCA analysis of the top 20 PCs was performed (Figure 6D). The t-SNE algorithm was used to cluster cells, obtaining 21 clusters (Figure 7A). We showed the top ten marker genes for the 21 clusters (Supplementary Table S5). In the Sepsis group, monocyte clusters, T Cells, and NK cells were decreased, and the B cell subpopulation increased (Figure 7B). We extracted mainly monocytes, NK cells, T Cells and B cells from the sepsis single cell dataset. GO enrichment analysis was performed to obtain the pathways of the differential genes (Supplementary Table S6-9) (Supplementary Figure S3-6). We show differential genes by plotting volcanoes, where red dots represent upregulated genes, blue dots represent downregulated genes, and yellow dots represent the seven core genes (Figure 7C). We found the most significant differences in CFLAR, STOM and BATF. Similarly, we compared the expression levels of the seven core genes in normal subjects and sepsis patients. In both groups, STOM, CASP4 and CFLAR were expressed at higher levels, while the remaining four genes were expressed at lower levels (Figure 7D).

A heatmap showing the proportion of core genes expressed in immune cells is then displayed. A high expression level of STOM, BATF, CASP4, and CLFAR was found in all samples in all 6 cell types (Figure 8A), and CLFAR and STOM were expressed at a high level at the gene expression level as well. STOM was highly expressed in platelet subpopulations in the normal and Sepsis groups (Figure 8B), while CLFAR expression was higher in the monocyte subpopulation, neutrophil subpopulation, and NK cell subpopulations in the Sepsis group compared to the normal group. We found that these results were generally consistent with what we found in Figure 6D in our previous analysis. Moreover, performing ssGSEA metabolic pathway analysis, we found that the glucose metabolism scores of monocytes and NK cells were different in normal and Sepsis (Figures 8C, D). Sepsis patients had higher glucose metabolism scores on monocytes and NK cells than normal patients, and CLFAR appears to be involved in this pathway (Figure 8E).