Three microarray datasets in this study were obtained from the Gene Expression Omnibus (GEO) database using the “GEOquery” package (version 4.2.1) (15). Among them, GSE9960 and GSE28750 both used the GPL570 sequencing platform (16). GSE9960 contained 70 peripheral blood mononuclear cell samples, including 54 sepsis samples as the case samples and 16 normal controls as the control samples. GSE28750 contained 30 blood samples, of which 10 sepsis samples were the case samples and the remaining were normal controls. The scRNA-seq dataset based on the GPL24676 sequencing platform, GSE167363 contained 12 peripheral blood mononuclear cells samples, including 10 sepsis case samples and 2 normal control samples (17). The 2086 telomere-related genes (TRGs) were extracted from TelNet database (http://www.cancertelsys.org/telnet) (18). The analysis process of this study is shown in Figure 1.

The PCA in GSE9960 dataset was performed through the “scatterplot3d” package (v 0.3-42) (19). The differentially expressed genes (DEGs) between sepsis and normal samples in GSE9960 dataset were identified using the “limma” package (v 3.54.1) (20). The screening criteria were P < 0.05 & |log₂Fold Change (FC)| > 0.5. The top 10 genes up-regulated and down-regulated in the sepsis samples based on the |log₂FC| values were displayed by volcano plot and heat plot which made use of “ggplot2” package (v 3.4.1) (21) and “ComplexHeatmap” package (v 2.14.0) (22), respectively.

In our study, the relative distribution ratios of 22 immune cells in each sample in GSE9960 were analyzed by using the cell identification by estimating relative subsets of RNA transcripts (CIBERSORT) algorithm (v 0.1.0) (23), and immune cells with a proportion of 0 over 30% samples were removed in the subsequent analysis. The differences of immune cells infiltration between the case and control group were completed through the Wilcoxon rank sum test (P < 0.05).Module genes highly correlated with immune cells were selected using the “WGCNA” package (v 1.71) (24). Sample cluster analysis used the hierarchical clustering (HCLUST) function to identify and eliminate outliers. Then, we set R² = 0.85 to screen for soft thresholds (β). Topology overlap and adjacency matrices were established based on the gene expression data. Afterwards, the minimum gene count per module was set to 100 with the altitude of models being set to 0.4, and the modules were merged according to the standards of the dynamic tree cutting algorithm. Then, differential immune cell scores were used as the phenotype to calculate the correlations between key module and phenotype. Finally, the modules with the best comprehensive score from the modules which significantly correlated with phenotype were obtained as the final module. The genes in final modules significantly correlated with phenotype were selected as the immune cell-related genes (ICRGs) (P < 0.05).

The candidate genes in our study were obtained by overlapping the ICRGs, DEGs and TRGs, which made use of the “ggvenn” package (v 1.7.3) (23).

The functions of candidate genes were explored by the kyoto encyclopedia of genes and genomes (KEGG) and gene ontology (GO) analyses with the “clusterProfiler” package (v 4.2.2) (25) (P < 0.05). Furthermore, the protein level interactions of candidate genes were analyzed by PPI network using data from the search tool for the recurring instances of neighboring genes (STRING) database (confidence > 0.4) (26).

In order to identify feature genes associated with sepsis, the leave-one-out cross-validation (LOOCV) framework was utilized in the GSE9960 and GSE28750 datasets that integrated ten different machine learning algorithms into 101 algorithm combinations. The input data for machine learning algorithms was the candidate genes, and the ten machine learning algorithms comprised RSF and Enet by the “glmnet” package (v 4.1.4) (27), least absolute shrinkage and selection operator (LASSO) by the “glmnet” package (v 4.1.4) (27), Ridge by the “glmnet” package (v 4.1.4) (27), stepwise Cox by the “caret” package (v 6.0-93) (28), xBoost by the “xgboost” package (v 1.7.3.1) (29), plsRcox by the (v 3.2.2), SuperPC by the “superpc” package (v 1.12) (30), generalized linear mixed model (GBRM) by the “gbm” package (v 2.1.8.1) (31), and survival-SVM by the “e1071” package (v 1.7-13) (32). Subsequently, the receiver operating characteristic (ROC) curve drawn using the “pROC” package (v 1.18.0) (33) was used to validate the predictive performance of the 101-machine learning algorithm combinations when the minimum and the maximum number of genes in the model was from 2 to 10. The combination model with the highest area under the curve (AUC) value of ROC curve was considered as the optimal model. The genes of the optimal model were considered as the candidate biomarkers for subsequent analysis.

In this study, ROC curve and expression level analysis were also performed on the GSE9960 and GSE28750 datasets to identify biomarkers with sepsis (AUC > 0.7). The expression level analysis was used by Wilcoxon test, which required that the expression trends of genes in two datasets were consistent (P < 0.05).

In this study, the genotype-tissue expression (GTEx) database was employed mainly for querying the expression levels of biomarkers in 28 normal tissues and organs of the human body. Besides, correlations among biomarkers were carried out using Spearman’s correlation test with the “psych” package (v 2.2.9) (34) in the GSE9960 datasets.

To investigate the specific role of biomarkers in the diagnosis of sepsis, the nomogram based on biomarkers was constructed using the “regplot” package (v 1.1) (35) in the GSE9960 datasets. The accuracy of prediction with nomogram was determined by the calibration curve, which was drawn by the “rms” package (v 6.5.0) (36). Furthermore, the predicting value of the nomogram was also appraised by ROC curve.

In order to explore the biological pathways and functions with biomarkers, gsea was executed by “clusterProfiler” package in the GSE9960 datasets (37). Primarily, the Spearman correlation coefficients between biomarkers and all genes were calculated in the disease samples of the GSE9960 dataset. The GSEA of biomarkers was conducted according to the sequencing list of genes which corresponded to biomarkers (P.adjust < 0.05, |normalized enrichment score (NES) |> 1). In this step of analysis, the gene set from the “org.Hs.eg.dbR” package (v 3.1.0) (38) was employed as the background gene set.

In order to cluster the sepsis samples into different clusters in the GSE9960 dataset, the k-means algorithm with 1000 iterations was executed by “ConsensusClusterPlus” package. In order to determine the expression of biomarkers and the differences in immune cell infiltration different subtypes on sepsis, the expression of biomarkers and immune infiltration analysis was performed between the different clusters according to the previously referenced method (P < 0.05).

LncRNAs can control the expression of mRNAs by binding to shared miRNAs. So, in this study, the miRTarBase v9.0 database and TarBase v9.0 database were manipulated based on the NetworkAnalyst platform to forecast the miRNAs which could target biomarkers. Then, the final miRNAs were obtained by crossing over the miRNAs from the two databases. Finally, the miRNet database was used to predict the lncRNAs which could target the final miRNAs, and the lncRNA-miRNA-mRNA network was visualized by Cytoscape software (v 3.9.1) (39).

In addition, this study also employed ChEA3 database to forecast the transcription factors (TFs) which could target biomarkers. The biological function (P < 0.05) and the distribution in the tissues of TFs also were obtained from the ChEA3 database.

To search for potential therapeutic drugs related to biomarkers for sepsis, the “enrichR” package (v 3.2) (39) was used based on the Drug Signatures Database (DSigDB) database to predict genes-drug interactions (P < 0.05). Afterwards, in order to evaluate the binding ability between drugs and biomarkers, the drugs which had the highest significance with biomarkers were selected to perform molecular docking. The three-dimensional structure with proteins of biomarkers was retrieved from the Protein Data Bank (PDB) database (https://www.rcsb.org/) and the three-dimensional structure with molecular ligands of key active ingredients was retrieved from the PubChem database. Finally, molecular ligands and proteins were docked using online website.

The scRNA-seq analysis was conducted using “Seurat” package (v 5.0.1) (40). To ensure the accuracy of single-cell data, all samples in the GSE167363 dataset were dealt with the PercentageEigenSet function of the “Seurat” package. The data processing conditions were as follows: the number of genes in cells ranges from 300 to 10000, the expression level of genes in cells between 300-2000, the genes expressing in at least three cells, and the proportion of mitochondrial genes less than 15% in cells. Then, in light of the GSE167363 dataset, data were normalized by the “NormalizeData” function in the “Seurat” package (v 5.0.1), and highly variable genes (HVGs) were selected by the “FindVariableFeatures” function. Next, the “ScaleData” function in the “Seurat” package (v 5.0.1) was applied to scale data before principal components analysis (PCA). Subsequently, the “JackStraw” function within the “Seurat” package (v 5.0.1) was applied to execute PCA on HVGs. The “ElbowPlot” function within the “Seurat” package (v 5.0.1) was thereafter applied to draw a scree plot of the top 30 principal components (PCs), aiming to identify PCs that notably contributed to variation for subsequent analysis (P < 0.05). Afterward, cell cluster analysis was conducted on cells after dimensionality reduction utilizing “FindNeighbors” and “FindClusters” functions (resolution = 0.1, dimension = 30). After that, the cells were clustered by the uniform tSNE clustering method (41). Marker genes for cell annotation in this study were obtained from the literature (15).

To determine the key cell types of sepsis in the GSE167363 dataset, the proportion of different types of cells in different samples, the differential infiltration of cells, and expression of biomarkers in different cell types were all considered. So, the cells that satisfied the following criteria simultaneously as key cells: (1). high proportion of cells in the sample; (2). cells with different infiltration ratios between the disease and control; (3). cells with different expression of biomarkers (P < 0.05). Subsequently, the secondary clustering analysis on key cells was performed which referred to the previous method, and the pseudo time analysis was proceeded by the “Monocle” package (v 2.30.0) (42) to study the developmental trajectory and differentiation directions of key cell types. Eventually, the expression of biomarkers during cell differentiation was observed.

The assessment of biomarkers expression was conducted on clinical tissue samples using RT-qPCR. A total of 5 pairs of blood samples were obtained from Wuhan Central Hospital Affiliated to Huazhong University of Science and Technology, including 5 sepsis and 5 control. All participants needed to sign and fill the informed consent form, and the ethical approval agency was Ethics Committee of The Central Hospital of Wuhan (No. WHZXKYL2024-164). Firstly, the total RNA of 5 pairs of tissue samples was derived by TRizol reagent (Ambion, U.S.A). The RNA concentrations were computered by NanoPhotometer N50. Secondly, mRNA was reversely transcribed into complementary DNA (cDNA) utilizing SureScript-First-strand-cDNA-synthesis-CREB5B test kit (Servicebio, Wuhan, China). Finally, the RT-qPCR was conducted. The expression levels of biomarkers between sepsis and control samples were calculated by 2^-ΔΔCt. The internal reference gene was glyceraldehyde-3-phosphate dehydrogenase (GAPDH), which was employed to normalize the results. The results were calculated by GraphPad Prism 5. Detailed information of primers and machine testing conditions was listed in Supplementary Table S1.

Bioinformatics analyses were performed utilizing the R programming language (v 4.2.2). Wilcoxon rank sum test was used to compare the differences between two groups. P < 0.05 was considered statistically significant.

The results of PCA showed that the case and control samples in the GSE9960 dataset were relatively separated, but there were no outlier samples, indicating good stability of the group samples (Supplementary Figure S1). So, the DEGs were obtained between all case and control samples in this dataset. A total of 589 DEGs were identified, of which 342 were up-regulated and 247 were down-regulated in the case group (Figures 2A, B). In the immune infiltration analysis, 15 types of immune cells were remained after removing cells that did not conform to the criteria (Supplementary Figure S2). Then, a significant difference was observed in the infiltration of 10 immune cells between the case and the control group, including memory B cell, plasma B cell, Macrophage M0, Monocyte, Neutrophil, resting natural killer (NK) cell, CD+ resting memory T-cell, naive CD4+ T cell naive, CD8+ T cell CD8+, and gammadelta T cell. The 10 immune cells were included in subsequent analysis (Figure 2C).

In WGCNA, no outlier samples were detected and subsequent analysis was based on all samples (Figure 2D). Afterwards, the soft threshold was 30 when R² = 0.85, and a total of 7 co-expression modules were obtained (Figures 2E, F). The correlation analysis indicated that the MEyellow module exhibited the strongest correlation with the immune cell score, and 1,886 ICRGs were identified (Figure 2G). Ultimately, 32 candidate genes were obtained by taking the intersection of ICRGs, DEGs, and TRGs (Figure 2H).

The GO enrichment analysis of the candidate genes revealed that there were 440 significant functions, including nuclear division, organelle fission, spindle, chromosomal region, protein serine kinase activity, protein serine/threonine kinase activity, etc. (Figure 3A) (Supplementary Table S2). Meanwhile, a total of 15 significant functions were enriched in the KEGG analysis which included oocyte meiosis, progesterone-mediated maturation of oocytes, cellularsenescence, cell cycle, vral carcinogenesis, etc. (Figure 3B) (Supplementary Table S3). The PPI network was established based on the candidate genes, 19 genes had interactive network relationship. It was worth noting that CDK1 exhibited interactions with the most genes (Figure 3C).

A total of 36 models were obtained when machine learning combination models were selected, among which the “RF+Lasso” model was considered the optimal model due to its highest value of AUC (Figures 4A–C). At the same time, the 7 genes such as MYO10, TDRD9, SULT1B1, MKI67, CREB5, BASP1, and CKAP4 in the “RF+Lasso” model were placed in subsequent analysis. The ROC curve analysis of the7 genes indicated that the AUC values of MYO10, SULT1B1, MKI67, and CREB5 were greater than 0.7 in the two datasets (Figures 4D, E). These 4 genes were differentially expressed between groups in both datasets, and all were up-regulated in the case group (Figures 4F, G). Therefore, MYO10, SULT1B1, MKI67, and CREB5 were selected as the biomarkers for this study. At that time, all 4 biomarkers showed a positive correlation, and the highest positive correlation occurred between SULT1B1 and CREB5 (cor = 0.823182574, p = 2.24×10^-18) (Figure 4H). Otherwise, the GTEx database was occupied to observe the distribution of biomarkers in human tissues and organs (Figure 4I). The results indicated that CREB5 and SULT1B1 had the highest expression level in whole blood tissue, and MKI67 had the highest expression level in cell-Cultured fibroblasts. Unusually, the levels of expression with MYO10 in all tissues were similar.

The nomogram consisted of “points” and “total points”, with the points representing the points of biomarkers and the latter representing the total points of all biomarkers. The sum of gene points indicated that the higher the total score, the higher likelihood of sepsis in this sample. As shown in the Figure 5A, the point of MYO10 was the highest indicated that MYO10 had the greatest impact on the predictive value of sepsis and the result was significant (P < 0.05). In addition, the P-value of the calibration curve was 0.111, indicating the accuracy of prediction with the nomogram was quite well (Figure 5B). Similarly, the AUC in ROC of the model was 0.872, also indicating good prediction performance of the model (Figure 5C).

The GSEA revealed the enrichment pathway of the biomarkers. The results showed that the four biomarkers were all enriched in ribosome, and 3 biomarkers (MYO10, SULT1B1, CREB5) were enriched in Fc gamma R-mediated phagocytosis. The common enrichment pathway of CREB5 and SULT1B1 was the chemokine signaling pathway. The co-enrichment pathway of CREB5 and MYO10 was the regulation of actin cytoskeleton (Figures 6A–D).

As shown in the Figure 7A, the best cluster stability was achieved when the number of clusters (K) equaled 2. Consequently, the sepsis samples in the GSE28750 dataset were divided into two subgroups, including cluster1 and cluster2. To obtain the expression trends of biomarkers between cluster1 and cluster2, the Wilcoxon test was used. The results indicated that the expression of MKI67, SULT1B1 and CREB5 were significantly different between two clusters. Specifically, MKI67 was up-regulated in cluster2, while CREB5 and SULT1B1 were up-regulated in cluster1 (Figure 7B). To investigate the different proportion of immune cell infiltration between two clusters, 22 immune cells were included based on the CIBERSORT algorithm (Supplementary Figure S3). But only 16 immune cells were included in differential analysis, and only two cells stypes (B cell plasma b and Neutrophil) showed significant differences between clusters (P < 0.05) (Figure 7C).

By using the NetworkAnalyst platform, a total of 30 miRNAs were identified for four biomarkers. Among them, 3 miRNAs were predicted with CREB5, 17 miRNAs were predicted with MKI67, 8 miRNAs were predicted with MYO10, and 2 miRNAs were predicted with SULT1B1. Specifically, the miRNA which targeted MYO10 and SULT1B1 was hsa-mir-124-3p. while that targeted MYO10 and MKI67 was hsa-mir-218-5p. Next, 1,055 lncRNA-miRNA interactions were obtained by the miRNet database, and the top 3 lncRNAs related to each miRNA such as which included MKI67-hsa-mir-484-SNHG12, MKI67-hsa-mir-218-5p-SLFNL1-AS, and MYO10-hsa-mir-218-5p-ADAMTSL4-AS1 were selected to construct a network diagram. The lncRNAs-miRNAs-biomarkers network were exhibited in the Figure 8A, Besides, a total of 703 TFs were identified for four biomarkers. We selected the top 30 TFs which had higher correlation with the biomarkers for the TF-biomarkers network. The network was visualized by Cytoscape software. The number of TFs predicted with MKI67 was the largest (Figure 8B). The enrich function of TF mainly included anatomical structure morphogenesis, animal organ morphogenesis, and blastoderm segmentation. Finally, we also observed the expression of TFs in the tissue distribution, and the results displayed that the expression of TFs were observed in nerve, colon, uterus, skin, and blood vessel (Figure 8C).

A total of 80 drugs were chosen in the DSigDB database. Subsequently, we selected the top 10 drugs with higher significance to construct drug-biomarkers network by Cytoscape software. In the drug-biomarkers network, it could be observed that Methaneseleninie acid was targeted to the MYO10 and MKI67 while MS-275 was targeted to the MYO10 and CREB5 (Figure 9A). Due to its highest saliency with biomarkers and three-dimensional structure, MS-275 was selected to perform molecular docking analysis. MYO10 was selected as the gene also because of the three-dimensional structure for molecular docking analysis. It is generally believed that |total score| > 7.0 kcal/mol indicated a stronger binding activity. In our study, the total score was -9.4 kcal/mol, indicating that the biomarkers had preferable binding activity with the target protein (Table 1) (Figure 9B).

After undergoing quality control, 47,080 cells and 20,696 genes were retained for further analysis (Supplementary Figure S4). Afterwards, the top 3000 hypervariable genes were selected for PCA (Figure 10A). In this study, the top 30 principal components were selected for cluster analysis, and all cells were ultimately divided into 12 clusters (Figures 10B, C). A total of 6 cells (B cells, CD16+ and CD14+ monocytes, CD4+memory cells, CD8+ T cells, Megakaryocyte progenitors, NK cells) were annotated based on the expression of marker genes in different cell clusters (Figure 10D). Simultaneously, the expressions of marker genes were shown in the six cell types (Figure 10E). To further demonstrate the expression of marker genes in each cell cluster, we generated a heatmap of the top 5 marker genes with the highest expression levels in each cluster (Supplementary Figure S5), aiming to more intuitively reflect the gene expression characteristics of different cell clusters. The cell proportion bar stack graph shows the high specificity of the marker gene, confirming the accuracy of the annotation (Figure 10F). To identify the key cell type in this study, we analyzed the expression levels of biomarkers in different cells. The results showed that compared with other immune cell types, all four biomarkers exhibited high expression in CD16+ and CD14+ monocytes (Figure 10G). We further compared the infiltration differences of each cell type between the sepsis group and the control group using the Wilcoxon test (P < 0.05), and found that there were significant differences in CD4+ memory cells, B cells, CD16+ and CD14+ monocytes, and CD8+ T cells between the two groups (Figure 10H). By synthesizing the evidence from both the above-mentioned biomarker expression profiles and cell infiltration differences, we finally identified CD16+ and CD14+ monocytes as the key cells. After secondary clustering, key cells were divided into seven clusters of key cells (Supplementary Figure S6). In pseudo-time analysis, a total of 5 underwent distinct states of CD16+ and CD14+ monocytes were discovered, and the expression of the cell was reduced over time in the disease group (Figure 11A).

During the process of cell differentiation, MKI67 and MYO10 showed no significant changes, while the expression level of CREB5 first decreased and then increased, mostly distributed in the state 1 stage, which was the early stage of differentiation. The expression levels of SULT1B1 first decreased and then stabilized, mostly distributed in the state 1 stage. The results indicated that the development of diseases might be closely related to CREB5 and SULT1B1 (Figure 11B).

As shown in the Figure 12, The RT-qPCR results showed the expression of MYO10, SULT1B1, MKI67, and CREB5 had significant differences between controls and sepsis samples (P < 0.05). The expression of 4 biomarkers were all upregulated in the sepsis group which were consistent in the results with the bioinformatics analysis results, indicating that preliminary results were reliable in our study.