Identification of key genes as potential diagnostic biomarkers in sepsis by bioinformatics analysis

Dataset acquisition and combination

Sepsis-related gene expression profile datasets were obtained from the Gene Expression Omnibus (GEO) database (Clough & Barrett, 2016). The GSE28750 and GSE74224 datasets were used for exploratory analyses. The GSE28750 dataset consisted of 10 sepsis cases and 11 critically ill controls. The GSE74224 dataset included 74 sepsis cases and 31 critical illness controls. Table S1 showed the overview of sepsis datasets used in this study. Both datasets compared ICU sepsis patients and post-surgical patients with infection-negative systemic inflammation (McHugh et al., 2015; Sutherland et al., 2011). The “removeBatchEffect” function of the “limma” R package was used to remove the batch effect of the two datasets, thus merging the two datasets into one integrated gene expression dataset (Ritchie et al., 2015). All subsequent analyses were based on the merged dataset.

Weighted gene co-expression network analysis

Weighted gene co-expression network analysis (WGCNA) was performed to identify gene modules highly associated with sepsis by utilizing “WGCNA” R package (Langfelder & Horvath, 2008). Firstly, the “pickSoftThreshold” function determined the ideal weighted parameters for adjacent functions, which then served as soft thresholds in the construction of the network. Following this, a weighted adjacency matrix was established. To categorize gene modules, hierarchical clustering was employed, leveraging a topological overlap matrix (TOM) and employing a dissimilarity index (1-TOM). The subsequent step involved evaluating the association of each gene module with the specific disease, identifying the module with the strongest correlation as the key module. Genes within this key module were then earmarked for further investigation.

Enrichment analyses of key module genes

The key module genes identified by WGCNA were subjected to Gene Ontology Biological Process (GO_BP), Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome enrichment analyses via the Database for Annotation, Visualization and Integrated Discovery (DAVID; https://david.ncifcrf.gov/) (Sherman et al., 2022). The results of the enrichment analysis were visualized by plotting dot plots using the “ggplot2” R package, exhibiting the top 10 most significantly enriched terms in each category.

Differential expression analysis

Differential expression analysis was performed using the “limma” R package (Ritchie et al., 2015). Subsequently, differential expression analysis was performed using “lmFit”, “eBayes”, and “topTable” functions. Genes with adjusted P-value (false discovery rate (FDR)) < 0.05 and —log₂ fold change— (—logFC—) > 0.5 were identified as differentially expressed genes (DEGs). Finally, the DEGs were visualized using the “ggplot2” R package and the “pheatmap” R package to respectively generate volcano maps and heat maps.

Construction of protein-protein interaction networks for DEGs

The protein-protein interactions (PPI) network of DEGs was constructed through the STRING online platform (https://david.ncifcrf.gov/) (Szklarczyk et al., 2023). Subsequently, the PPI network was imported into Cytoscape software and screened for hub genes based on the cytoHubba plugin (Chin et al., 2014). A total of four diverse cytoHubba algorithms, namely Degree, EPC, MCC, and MNC were used for screening hub genes.

Identification of sepsis-related key genes

Using the “UpSetR” R-package and the “VennDiagram” R-package, the intersecting genes of both the genes in the key modules identified by WGCNA and the hub genes under the four algorithms were obtained as the sepsis-related key genes.

Validation of key genes

Validation of key genes was performed in an independent dataset, GSE131761 (containing 81 sepsis cases and 33 critical controls; Table S1). Specifically, the ability of each key gene to distinguish among sepsis cases and critical controls was first identified in the GSE131761 dataset using the “pROC” R package with a receiver operating characteristic curve (ROC). Subsequently, the genes in the WGCNA key module were defined as the gene set for sepsis signature and the sepsis signature level of each sample in the GSE131761 dataset was quantified using the “ssgsea” algorithm of the “GSVA” R package. Finally, based on the results of the normality test, the differences in sepsis signature levels and key gene levels between sepsis cases and critical controls were assessed using the t-test or Wilcoxon test.

Immune cell analysis

The relative proportion of immune cells for each sample in the combined dataset was assessed using the “CIBERSORT” algorithm (Chen et al., 2018). Subsequently, based on the results of the normality test, a t-test or Wilcoxon test was used to compare the differences of each immune cell between sepsis cases and critical controls. Subsequently, based on the results of the normality test, Pearson or Spearman correlation analyses were performed to assess the correlation of each key gene with the immune cells. The results of the correlation analysis were visualized using the heatmap.

Validation in sepsis patients

To further corroborate the effectiveness of the key genes, we recruited 22 sepsis patients and 10 non-sepsis critically ill controls from the intensive care units (ICUs) in our hospital between September 2020 and December 2021. All the participants were 18 years or older and had an ASA classification of III-IV. Sepsis patients were included when they met the diagnostic criteria for sepsis (Singer et al., 2016) and had a clear focus of infection in the urinary (obstructive urinary tract disease), gastrointestinal (gastrointestinal tract perforation) or biliary tracts (acute purulent obstructive cholangitis). Non-sepsis critically ill controls were patients with multiple traumatic injuries of similar age to those with sepsis and without co-infection. Patients were excluded if they received antibiotics or hormones before ICU admission; were currently receiving renal replacement therapy; had any cancers, autoimmune diseases, cirrhosis, uremia or had recent history of blood transfusion. Comparison of basic information for patients with sepsis and non-sepsis critical illness was shown in Table S2.

All blood samples were collected within 24 hours of ICU admission and prior to emergency surgery, after which PBMCs were isolated. RNA was then extracted from the PBMCs with TRIzol reagent (12183-555; Invitrogen) according to the manufacturer’s instructions. cDNA was synthesized using the SuperScript™ III First-Strand Synthesis SuperMix (11752-050; Invitrogen) for quantitative real-time polymerase chain reaction (qRT-PCR). The qRT-PCR was performed utilizing the HiScript II Q RT SuperMix for qPCR (Vazyme, Nanjing, China) in a LightCycler 480 PCR system (Hoffmann-La Roche Ltd., Shanghai, China). The relative mRNA expression levels were computed through the application of 2 - ΔΔCT method, with GAPDH serving as an internal reference. The primer pairs used in the experimental procedure are listed in Table S3. The data was presented as mean ± SD and analyzed via the Student’s t-test. Statistical significance was established when the p value was less than 0.05.

Our study was approved by the Institutional Review Board of the Third Xiangya Hospital (2020-S373) in strict compliance with the Declaration of Helsinki. Furthermore, written informed consents were procured from all patients or their legally designated representatives in our study.

Dataset combination and DEGs identification

The overall workflow of this study is illustrated in Fig. 1. Principal Component Analysis (PCA) indicated that samples from GSE28750 and GSE74224 were batch-effect corrected and formed a new combined dataset (Figs. 2A–2B). Subsequently, in the combined dataset, a total of 293 DEGs (including 173 upregulated DEGs and 120 downregulated DEGs) were identified between sepsis and critical control (Figs. 2C–2D).

WGCNA and identification of key module

First, we opted for a β value of 5, corresponding to a scale-free R² of 0.85, to establish scale-free networks (Fig. 3A). After merging similar modules, a total of 9 gene modules were identified in the combined dataset (Fig. 3B). Correlation analysis revealed that the greenyellow module was significantly correlated with sepsis (correlation coefficient: 0.7, P = 2e−19) (Fig. 3C). Genes within the greenyellow module will be used for subsequent analysis. Subsequently, correlation analysis between module membership and gene significance revealed that these genes were highly correlated with both module and phenotype (cor = 0.64, p = 2.3e−128) (Fig. 3D).

Enrichment analysis of greenyellow module genes

The greenyellow module genes were subjected to enrichment analysis through the DAVID platform, and the results were presented in dot plots (Fig. 4). Specifically, upon conducting GO BP enrichment analysis, it was observed that genes within the greenyellow module showed significant enrichments in processes such as apoptotic process, cell division, intracellular signal transduction, protein phosphorylation, and inflammatory response (Fig. 4A). As for the KEGG pathway enrichment analysis, these genes exhibited notable associations with pathways including MAPK signaling pathway, prion disease, diabetic cardiomyopathy, chemical carcinogenesis—reactive oxygen species, and cell cycle (Fig. 4B). Moreover, analysis of Reactome pathways highlighted significant enrichments encompassing immune system, innate immune system, cell cycle, and neutrophil degranulation among the genes within the greenyellow module (Fig. 4C).

Construction of PPIs for DEGs and identification of hub genes

Through the STRING platform, a PPI network of DEGs was constructed, containing 232 nodes and 821 edges (Fig. 5A). The PPI network included 143 up-regulated DEGs and 89 down-regulated DEGs (Table S4). In addition, based on four different algorithms, namely Degree, EPC, MCC, and MNC, top 20 hub genes in the PPI network were identified for subsequent analysis (Fig. 5B).

Identification and validation of sepsis key genes

Upset plot and Venn diagram showed that LTF, LCN2, ELANE, MPO, and CEACAM8 existed both in the greenyellow module and in the top 20 hub genes under the four algorithms (Fig. 6A). In the external validation dataset GSE131761, the ROC demonstrated that five key genes exhibited a promising ability to differentiate sepsis patients from critical controls (Fig. 6B). Using the “ssGSEA” algorithm, sepsis signature scores were found to be significantly higher in sepsis patients than in critical controls in the GSE131761 dataset, thus demonstrating the reliability of the greenyellow module identified by the WGCNA (Fig. 6C). In addition, the expression levels of LTF, LCN, ELANE, MPO, and CEACAM8 were all up-regulated in the peripheral blood of sepsis patients than critical controls in the GSE131761 dataset (Fig. 6D), which is consistent with the previous findings.

Subsequently, based on an external validation cohort from our hospital, qRT-PCR demonstrated significant up-regulation on LTF, LCN2, ELANE, MPO, and CEACAM8 in peripheral blood of sepsis patients versus controls (Fig. 7), thus further enhancing the reliability of the findings.

Immune cell analysis

Figure 8A shows the distribution of circulating immune cells in the combined dataset of 84 sepsis patients and 42 critical controls. Figure 8B shows that circulating monocytes were significantly upregulated in sepsis patients compared to critical controls, whereas mast cells resting and T cells CD8 were significantly downregulated compared to critical controls. Correlations between immune cells and five key genes were assessed in patients with sepsis, revealing that plasma cells, macrophages M0, monocytes, T cells regulatory (Tregs), eosinophils and NK cells resting were simultaneously and significantly associated with more than two key genes (Fig. 8C).