All GEO datasets were downloaded from Gene Expression Omnibus (GEO, www.ncbi.nlm.nih.gov/geo/) (Barrett et al., 2013). We selected 6 datasets (GSE154918, GSE54514, GSE9960, GSE69063, GSE25504, and GSE13904) related to sepsis for analysis. The whole gene expression profiles of peripheral blood were extracted for further analysis. The GEO datasets collected are exhibited in Table 1. The flowchart of the study was elucidated in Figure 1.

The GSE154918 dataset (GPL20301 platform) is composed of whole gene expression profiles of peripheral blood from 40 control and 24 sepsis samples, The GSE54514 dataset (GPL6947 platform) was composed of whole gene expression profiles of peripheral blood from 36 non-sepsis healthy control subjects and 127 sepsis patients, The GSE9960 dataset (GPL570 platform) was composed of whole gene expression profiles of peripheral blood from 16 control and 54 sepsis samples, The GSE69063 dataset (GPL20301 platform) was composed of whole gene expression profiles of peripheral blood from 33 control and 57 sepsis samples, The GSE25504 dataset (GPL570 platform) was composed of whole gene expression profiles of peripheral blood from 37 control and 26 sepsis samples, The GSE13904 dataset (GPL570 platform) was composed of whole gene expression profiles of peripheral blood from 18 control and 52 sepsis samples.

R software’s “limma” package was applied for identifying the DEGs between the sepsis and control samples in the GSE154918 and GSE25504 datasets, respectively (Ritchie et al., 2015). DEGs with |log2 fold change (FC)| > 0.5 and p < 0.05 were defined as statistically significant. Volcano plots and heatmaps of the identified DEGs were visualized using the “ggplot2” and “heatmap” R packages.

R software’s “WGCNA” package was used to construct the co-expression network of the GSE154918 and GSE25504 datasets, respectively (Langfelder and Horvath, 2008). In brief, we explored the association among different pairs of genes and weighted them based on the expression levels of related genes in control and sepsis samples. Afterward, we transformed the adjacency matrix into a topological overlap matrix (TOM) to further verify the gene connectivity in the network. Next, the hierarchical clustering method was conducted to determine remarkably interconnected clusters (modules) according to gene connectivity and covariance coefficients. We selected the best power value and established the correlated modules. Finally, we performed Pearson correlation analysis and defined the highest three gene modules as statistically significant based on the correlation coefficient.

The GSE9960, GSE13904, and GSE54514 datasets were log2 transformed. The robust multi-chip average (RMA) method was utilized to normalize gene expression, and the batch effect between different platforms was eliminated using the ComBat method based on the “sva” R package. Principal component analysis (PCA) was applied for evaluating the aggregation between different samples. Afterward, R software’s “ConsensusClusterPlus” package was carried out to perform unsupervised consensus analysis (Wilkerson and Hayes, 2010) among the 233 standardized sepsis patients. The number of re-samplings settled at 1,000, with each re-sampling containing 80% of the samples. The maximum number of clusters was set to 10, and the optimal k value was determined by the cumulative distribution function (CDF) index and the consensus matrix. Finally, t-SNE was carried out to validate the subtype assignments based on the gene expression profiles of the above sepsis patients.

The dataset was normalized and summarized using robust multi-chip average (RMA) implemented in the R package affy, and batch effects were corrected.

R software’s “GSVA” and “GSEABase” packages were utilized to validate the performance of gene sets (c2. cp.kegg.v7.4. symbols) among sepsis patients with different molecular subtypes, thus identifying the enriched pathways in each subtype. Gene sets with adjusted p-value < 0.05 were defined as significantly enriched gene sets.

The DEGs were screened using DEGs analysis and the WGCNA method in GSE154918 and GSE25504 datasets, and the specific genes were also identified in sepsis patients with distinct molecular subtypes. Eventually, a total of 48 core genes were identified by intersecting all the results. The “ggVennDiagram” package was used in generating the Venn diagrams of co-DEGs.

Spearman correlation analysis was performed to determine the correlation between core genes based on the gene expression profiles. The heatmap of the correlation coefficient among these hub genes was visualized using the “corrplot” R package, and the gene relationship network diagram with a correlation coefficient >0.9 was constructed using the “igraph” R package.

The Database for Annotation, Visualization, and Integrated Discovery (DAVID, https://david.ncifcrf.gov/summary.jsp), an online tool, was utilized to perform the GO enrichment analysis of co-DEGs (Huang da et al., 2009). Statistically significant GO terms (BP) with FDR <0.05 were screened and the results were visualized using the “GOplot” R package.

KEGG and Reactome enrichment analysis of these co-DEGs were conducted using the “clusterProfile” and “org.Hs.eg.db” R packages (Yu et al., 2012). The statistically significant enrichment pathways with adjusted p-value < 0.05 were defined and visualized using a bubble plot.

The STRING database (https://strin g-db.org/) was applied for constructing the PPI network of co-DEGs, and the core genes with a combined score of more than 0.4 were screened (Szklarczyk et al., 2015). The protein–protein interaction information was visualized utilizing the Cytoscape software (version 3.8.2).

LASSO is a regularization method with strong predictability, that is better than regression analysis when examining high-dimensional data (Bader and Hogue, 2003). A total of 303 samples including 70 health and 233 sepsis patients from GSE9960, GSE13904, and GSE54514 datasets were applied for LASSO regression analysis. As a training set, 70% of samples (N = 217, 58 controls, and 169 sepsis samples) were randomly selected. As a validation set, 30% of samples (N = 86, 22 controls, and 54 sepsis samples) were selected. The LASSO regression analysis was performed based on the expression profiles of co-DEGs using the “glmnet” R package. The predictive model constructed in the training set was verified in the testing set. In addition, receiver operating characteristic (ROC) curves were drawn in the validation set and train set using the “pROC” R package to assess the performance of the constructive model. AUC more than 0.75 was defined as a model with a high diagnostic value.

GSE154918 and GSE69063 datasets were selected as the validation set to verify the diagnostic value of these genes identified by LASSO in sepsis. Genes with AUC >0.9 in both the GSE154918 and GSE69063 were defined as the ultimate hub genes. Subsequently, we verified the expression of these hub genes in GSE154918 and GSE69063 datasets and in distinct molecular subtypes of sepsis.

The flowchart of our study is shown in Figure 1. To identify sepsis-related genes, we firstly analyzed the DEGs between sepsis and control samples in the GSE154918 dataset. A total of 1,671 up-regulated and 1,623 down-regulated genes were determined using the DEGs method (Figures 2A,B, Supplementary Figure S1A). Subsequently, we applied the WGCNA method to study the co-expression network in the GSE154918 dataset. When the best soft threshold power settled at 6, the mean connectivity was more effective (Figures 2C,D). A total of 13 co-expressed gene modules were created. among which the blue (4,036 genes, R = 0.87, p < 0.05), turquoise (4,618 genes, R = 0.83, p < 0.05), and yellow modules (1719 genes, R = 0.8, p < 0.05) had the most correlation with sepsis patients (Figure 2E, Supplementary Figure S1B ). In addition, the blue (R = 0.85, p < 0.05), turquoise (R = 0.95, p < 0.05), and yellow (R = 0.92, p < 0.05) modules were significantly correlated with module-related genes (Supplementary Figures S1C–H). Meanwhile, we identified 945 up-regulated and 793 down-regulated genes using the DEGs method in the GSE25504 dataset (Figures 3A,B, Supplementary Figure S2A). Based on the WGCNA method, a relatively high mean connectively was maintained by setting the best soft threshold power to 8 (Figures 3C,D). We clustered a total of 14 co-expressed gene modules, among which the blue (2,174 genes, R = 0.79, p < 0.05), brown (1717 genes, R = -0.83, p < 0.05), and turquoise (2,478 genes, R = 0.8, p < 0.05) modules were closely correlated with sepsis patients (Figure 3E, Supplementary Figure S2B). Moreover, blue (R = 0.88 p < 0.05), brown (R = 0.93, p < 0.05), and turquoise (R = 0.9, p < 0.05) modules were also remarkably correlated with module-related genes (Supplementary Figures S2C–H).

The GSE9960, GSE13904, and GSE54514 datasets were mixed into a combined dataset containing a total of 233 sepsis samples. These three datasets exhibited obvious separation before batch correction (Figure 4A). While the batch effect among these datasets from different platforms had been successfully eradicated after batch correction (Figure 4B). Then, we performed molecular subtypes analysis in a combined dataset based on the expression of genes via the “ConsensusClusterPlus” R package. Clustering results suggested that the classification was most reliable and stable when k = 3 (Figures 4C–E, Supplementary Figure S3). Consistently, the t-SNE confirmed that only cluster1, cluster3, and cluster4 could be significantly separated (Figure 4F). In total, 233 sepsis patients were classified into three subtypes, including cluster1 (n = 144), cluster 3 (n = 26), and cluster4 (n = 39) based on gene expression levels, which were chosen for subsequent analysis.

The DEGs among different molecular subtypes of sepsis patients were identified based on whole genome expression profiles. Eight hundred and sixty-seven (473 up-regulated and 394 down-regulated) DEGs were identified between cluster 1 and cluster 3 samples (Supplementary Figures S4A,B). The top 10 overexpressed and suppressed DEGs between cluster 1 and cluster 3 were exhibited in Table 2 . In addition, a total of 420 (237 up-regulated and 183 down-regulated) DEGs were screened after comparing cluster 1 with cluster 4 samples (Supplementary Figures S4C,D). The top 10 overexpressed and suppressed DEGs between cluster 1 and cluster 4 were exhibited in Table 3. Moreover, one thousand four hundred and forty-five (696 up-regulated and 749 down-regulated) DEGs were acquired after comparing cluster 3 with cluster 4 samples (Supplementary Figures S4E,F). The top 10 overexpressed and suppressed DEGs between cluster 3 and cluster 4 were exhibited in Table 4 . Overall, we identified 1803 specific genes among three subgroups.

GSVA analysis was performed among distinct molecular subtypes of sepsis patients based on the expression profiles of DEGs. In total, 103 differentially enriched gene pathways, comprising 59 pathways with activation and 44 pathways with suppression were screened when comparing cluster1 with cluster3 samples (Supplementary Table S1). The top 20 differentially enriched gene pathways between cluster1 and cluster3 are exhibited in a heatmap (Figure 5A). A total of 90 differentially enriched gene pathways, consisting of 42 pathways with activation and 48 pathways with suppression, were obtained when comparing cluster1 with cluster4 samples (Supplementary Table S2). The top 20 differentially enriched gene pathways between cluster1 and cluster4 are exhibited in a heatmap (Figure 5B). A total of 110 differentially enriched gene pathways were obtained, including 56 pathways with activation and 54 pathways with suppression after comparing cluster3 with cluster4 samples (Supplementary Table S3). The top 20 differentially enriched gene pathways between cluster3 and cluster4 are exhibited in a heatmap (Figure 5C). Finally, a total of 34 enriched gene pathways among clusters were identified and presented in a heatmap (Figure 5D).

The DEGs identified by the DEGs and WGCNA methods in GSE154918 and GSE54514 datasets were intersected with the specific genes among different molecular subtypes of sepsis patients. Eventually, a total of 40 co-DEGs were screened (Figure 6A). Spearman’s correlation analysis was performed to illuminate the correlation patterns among these 40 co-DEGs, and the area of the pie chart represents the exact value of correlation coefficients. We found a significant correlation among these co-DEGs. For example, GPR84 show strong antagonistic effects with GYG1 (coefficient = 0.968), FCGR1B (coefficient = 0.968), and ANKRD22 (coefficient = 0.912). Simultaneously, CARD11 could also present synergistic effects with GYG1 (coefficient = -0.908), GPR84 (coefficient = -0.91), and BLOC1S1 (coefficient = -0.918) (Figure 6B). These results suggested the regulatory balance among these co-DEGs. In addition, the gene relationship network indicated that genes with a correlation coefficient of more than 0.9 were as follows: IRAK3, IL18R1, ATP6V1CA, GYG1, ATP9A, NMT2, ITK, EIF4B, CCND2, CARD11, FBXO21, BLOC1S1, GPR84, LRG1, ANKRD22, FCGR1B (Figure 6C), indicating that the molecular subtypes of sepsis may be the results of multi-gene interactions. Subsequently, we assessed the enrichment pathways in which these co-DEGs involve. The GO enrichment analysis suggested that these co-DEGs were primarily enriched in immune-related biological functions and pathways, such as immune response, T cell receptor signaling pathway, T cell costimulation, positive regulation of NF-kappaB transcription factor activity 4, positive regulation of T cell proliferation (Table 5 and Figures 6D,E). KEGG pathway enrichment analysis indicated that co-DEGs chiefly participated in autoimmune and inflammation-related diseases including Graft−versus−host disease, Inflammatory bowel disease, and Rheumatoid arthritis. Additionally, Immune responses and hematopoietic-related signaling pathways such as T cell receptor signaling pathway and Hematopoietic cell lineage were also closely associated with these co-DEGs (Figure 6F). Consistently, Reactome enrichment analysis also revealed that these co-DEGs were mainly implicated in immunomodulation and signal transduction such as Translocation of ZAP-70 to immunological synapse, Generation of second messenger molecules, Costimulation by the CD28 family, TCR signaling, and Downstream TCR signaling (Figure 6G). Combining these results, it can be inferred that abnormal immune function may be the critical pathogenesis of sepsis. In addition, the PPI network of 40 co-DEGs, including 22 nodes and 27 edges, was analyzed by the STRING website and visualized by the Cytoscape software (Supplementary Figure S5).

To further screen out potential biomarkers for sepsis, we performed the LASSO regression analysis based on the expression profile of co-DEGs. We separated a total of 303 samples in the combined dataset (GSE9960, GSE13904, and GSE54514 datasets, 233 sepsis and 70 control samples) into a training set (70%) and a validation set (30%). Eventually, a total of 25 potential biomarkers with a non-zero coefficient were obtained (Figures 7A,B). The AUC of the 25-gene signature was 0.9051 in the training set and 0.7955 in the validation set (Figure 7C), which suggests that the 25-gene-based model may be able to correctly diagnose sepsis. Subsequently, we also utilized external datasets to ascertain the diagnostic value of these potential biomarkers. The value of the AUC of these potential biomarkers in the GSE154918 and GSE69063 datasets were exhibited in a scatter plot (Figure 7D). Genes with AUC >0.95 in both the GSE154918 and GSE69063 datasets were as follows: ANKRD22, GPR84, GYG1, BLOC1S1, CARD11, NOG, and LRG1 (Figure 7E). These results reveal that the 7 core genes identified have the most ability to differentiate sepsis from healthy controls.

To further demonstrate whether these 7 hub genes were worth using in clinical practice, we validated the expression of these core genes in different datasets. In GSE154918 and GSE69063 datasets, ANKRD22, GPR84, GYG1, BLOC1S1, and LRG1 were all highly expressed in sepsis patients, whereas NOG and CARD11 were dramatically decreased in sepsis patients (Figures 8A,B). Meanwhile, these 7 core genes were also subtype-specific biomarkers. Therefore, we speculated that these 7 hub genes have a high diagnostic ability in sepsis patients and are closely related to different subtypes of sepsis.