Identifying Potential Effective Diagnostic and Prognostic Biomarkers in Sepsis by Bioinformatics Analysis and Validation

Introduction

Sepsis is a common clinical emergency and critical illness, and is regarded as one of the life-threatening infectious diseases.¹ There are approximately 48.9 million septic patients and 11 million septic patients die worldwide each year, accounting for 19.7% of all the global deaths.²

From the perspective of immunology, the essence of sepsis is life-threatening multiple organ dysfunction caused by immune response and regulation disorder caused by severe infection.³ In other terms, sepsis progresses from excessive immune activation to widespread immunosuppression.⁴ What is more, the main cause of sepsis death is important organ dysfunction, which is characterized by impaired vascular barrier function, disordered immune regulation, excessive coagulation, impaired vascular regulation and ischemia, which further leads to multiple organ failure and death.³ Growing evidence supports immunosuppression as a major cause of sepsis-related deaths,⁵ and the main reason of immunosuppression is the death of immune cells, especially the apoptosis of T and B lymphocytes.^6–8

Sepsis is considered a medical emergency because it can rapidly develop into organ dysfunction and death even if active medical treatment is taken immediately.⁹ Although our understanding of the pathogenesis of sepsis has progressed, it is still a worldwide health challenge due to high morbidity and high mortality. In the last few years, many technologies have been developed to diagnose sepsis reliably and quickly, but none of them is widely used in clinical practice.^10,11

Bioinformatics analysis helps people to explore the molecular mechanisms underlying disease development and provides a new and effective method for identifying potential diagnostic biomarkers, assessing prognosis and redirecting treatment. Therefore, to better improve early diagnosis techniques and optimise treatment of sepsis, we use bioinformatics to identify diagnostic and prognostic relevant biomarkers and potential therapeutic targets for sepsis.

Materials and Methods

Microarray Data

We downloaded data from the Gene Expression Omnibus (GEO) database, the largest and most complete public gene expression data resource in the world today.¹² The array-based gene expression profile dataset GSE69528 included whole blood samples from 83 septic patients and 55 uninfected controls. The array-based gene expression profile dataset GSE13904 included whole blood samples from 52 septic children, 106 septic shock children and 18 normal children. The array-based gene expression profile dataset GSE25504 included whole blood samples from 25 sepsis-infected human neonates and 37 controls. The array-based gene expression profile dataset GSE26440 included whole blood samples from 98 septic shock children and 32 controls. The array-based gene expression profile dataset GSE28750 included whole blood samples from 10 septic patients and 20 controls. The array-based gene expression profile dataset GSE145227 included peripheral blood mononuclear cell samples from 10 septic children and 12 controls. The array-based gene expression profile dataset GSE54514 included whole blood samples from 26 septic survivors, 9 sepsis non-survivors and 18 controls.

Screening for Differentially Expressed Genes (DEGs) from Blood Microarray Datasets of Sepsis Patients

DEGs between sepsis patients and control individuals in GSE69528 were screened by the “limma” package¹³ and volcano plots of DEGs were plotted using the “ggplot2” package¹⁴ to show the differential expression. P < 0.05 and | log2FC |> 1 were considered statistically significant.

Screening of Immune-Related DEGs

From the Immunology Database and Analysis Portal (IMMPORT) database (https://www.immport.org/home), we downloaded totaling 2483 immune-related genes (IRGs). The overlapping genes between IRGs and DEG were identified as immune-related differential DEGs (IRDEGs) by Venn diagram analysis, and were analysed by R language and cluster analysis packages for subsequent Gene Ontology (GO), Kyoto Gene and Genome Encyclopedia (KEGG) pathway enrichment analysis.

Estimate the 22 Immune Cell Types in Sepsis and Screen the Genes Most Related to Immune Cells in the Sepsis Group

Cibersort R package was used for cellular immune infiltration analysis to obtain the immune cell content matrix.¹⁵ The Pearson algorithm of R software was applied to calculate the correlation between the cell content and gene expression in the sepsis group samples, and the genes with a correlation coefficient greater than 0.7 were further intersected with IRDEGs to obtain common genes.

Screening Key Genes According to Receiver Operating Characteristic Curve Analysis

We calculated the area under the receiver operating characteristic curve (AUC) for common genes by using the pROC package in R¹⁶ in the six datasets, and visualized the results by using ggplot2. We selected genes with an AUC greater than 0.8 in all six datasets as the key genes for subsequent analyses.

Expression of Key Genes

We used R software and the ComplexHeatmap package to calculate the expression difference of key genes between septic patients and controls in 6 datasets. Further, we used Spearman analysis and Circlize package in R to calculate the correlation between key genes and then visualized the results.

Correlation Analysis of Key Genes and Immune Infiltrating Cells

We performed Spearman correlation analysis of key genes and immune infiltrating cells using R software and the ggstatsplot package, and visualised the results using the ggplot2 package.

TF-miRNA-mRNA Interaction Analysis

We obtained miRNA-key gene interactions from TarBase¹⁷ and mirtarbase¹⁸ to identify regulatory miRNAs that affect key genes at the post-transcriptional level. Meanwhile, in order to determine the regulatory TF controlling miRNA, TF-miRNA interaction was obtained using Jaspar database.¹⁹ All of the above use NetworkAnalyst to analyze topology parameters (ie, degree and medium centrality).²⁰ We further used Cytoscape software²¹ to integrate the miRNA-target gene network and the TF-miRNA network to construct a TF-miRNA-target network to discover major regulatory genes.

Associations Between Key Genes and Diseases

The Comparative Toxicogenomics Database (CTD; http://ctdbase.org/) is a public database describing gene–disease interactions and providing information on gene–disease interactions.²² Ten key genes are input into the gene search box of CTD database one by one, and the diseases associated with each key gene can be obtained respectively. We therefore use these data to identify and analyse associations between key genes and diseases, and to understand the diseases in which key genes primarily act.

Functional Enrichment in Key Gene

The AmiGO database (v2.0; http://AmiGO.geneontology.org/AmiGO/) is a web application that is primarily used to query, browse, visualize and download gene ontologies and annotations.²³ We used this database to confirm the enrichment of key sepsis-related genes that have been identified in this study.

Predictive Value of Key Genes in Sepsis Death

Using the data of the dataset GSE54514, the selected 10 key genes, gender, age, Acute Physiology and Chronic Health Evaluation II (APACHE II) score, and neutral percentage were used to perform multivariate and univariate regression analysis on the prognosis of death. APACHE II is one of the ICU scoring systems, higher scores suggest an increased severity.

Patient Enrollment and Sample Collection

From March 23, 2022 to March 30, 2022, patients diagnosed with sepsis in the Department of critical medicine of the First Affiliated Hospital of Chongqing Medical University were included in the study. We collected samples at the first time of admission (generally within 1 hour after admission) of patients with sepsis before the use of antibiotics. The follow-up treatment was carried out according to the routine treatment scheme. This part of the study was approved by the ethics committee of the First Affiliated Hospital of Chongqing Medical University (ethic reference number: 2022-02). Written informed consent was obtained from all donors in accordance with the Declaration of Helsinki and under Institutional Review Board–approved protocols. The inclusion criteria were: (1) the patient met the diagnostic criteria of sepsis 3.0,²⁴ (2) the patient was aged between 18 and 90 years, and (3) the patient was hospitalized for more than 24 hours. The following patients were excluded: (1) using drugs that affect immunity, (2) having malignant tumors, and (3) pregnancy or lactation.

Real-Time PCR Analysis

Total RNA was extracted from whole blood using trlzogen reagent (Tiangen biochemical technology (Beijing) Co., Ltd). We use high-capacity cDNA reverse transcription kit (TaKaRa TB Green Premix Ex Taq II, RR820A) and first-strand cDNA synthesis kit (TaKaRa PrimeScript RT reagent Kit, RR0474) for RNA reverse transcription and first-strand DNA synthesis. Specific PCR primer pairs are listed in Table 1. GAPDH was used for internal control, and 2–ΔΔCT method was used to calculate relative mRNA expression level.

Table 1 Specific PCR Primer Pairs

Statistical Analysis

Cox regression analysis was conducted at both univariate and multivariate levels to analyze independent factors affecting prognosis by using R software and Survival package. Statistical analysis and visualisation were performed using R3.6.3 and ggplot2. Statistical significance was calculated by Student’s t-test, Mann–Whitney U-test, chi-square test or Fisher’s exact test (as appropriate). Statistical tests were bilateral tests, and P values less than 0.05 were considered statistically significant.

Results

Identification of DEGs

The flow chart of this study is shown in Figure 1. In the dataset GSE69528, 1572 DEGs were screened between sepsis patients and controls in total, of which 755 were up-regulated genes and 807 were down-regulated genes (Figure 2A). There are 164 overlapping genes (IRDEGs) between IRGs and DEGs. The Venn diagram is shown in Figure 2B.

Figure 1 Flow chart of this study.

Figure 2 Identification and functional analysis of differentially expressed genes (A) volcano map, (B) Venn of GSE69528 and IMMPORT database, (C) BP, CC, MF analysis, (D) KEGG analysis).

GO and KEGG Pathway Enrichment Analyses of IRDEGs

IRDEGs obtained in the previous step are used for GO and KEGG analysis. BP associated with GO terms were mainly enriched in T cell activation, neutrophil degranulation, neutrophil activation, lymphocyte differentiation. CC were primarily associated with vesicle lumen, external side of plasma membrane, cytoplasmic vesicle lumen, secretory granule lumen. MF were mainly involved in receptor-ligand activity, cytokine receptor activity, MHC protein complex binding (Figure 2C).

KEGG pathway enrichment analysis of IRDEGs was mainly enriched in T cell receptor signaling pathway, cytokine–cytokine receptor interaction, Th17 cell differentiation, natural killer cell mediated cytotoxicity, Th1 and Th2 cell differentiation and Primary immunodeficiency (Figure 2D).

Immune Infiltration Analysis

To estimate 22 immune cell types, we used the Cibersort R software package to reveal the infiltration of 22 immune cell subsets in sepsis by analyzing the GSE69528 dataset. The box plot of the immune cell difference showed that, compared with the control sample, innate immune cells (such as mast cells, macrophages and neutrophils) infiltrated more, while T cells (CD4, CD8), NK cells infiltrated less, as shown in Figure 3A. Figure 3B shows the patterns of 22 subpopulations of immune cells’ proportion in sepsis and control group. The heat map of immune cell infiltration between the septic group and the control group is shown in Figure 3C.

Figure 3 Immune infiltration analysis. (A) Proportional box plot of 22 kinds of infiltrating immune cells in the sample. (B) The proportion and distribution of infiltrating immune cells in different samples. (C) Proportional heat map of 22 kinds of infiltrating immune cells in the sample.

To screen the most relevant genes of immune cells in sepsis group, the Pearson algorithm of R software was used to calculate the correlation between the content of immune cells and gene expression in the sepsis group samples. A total of 441 genes with a correlation coefficient greater than 0.7 were obtained, and a total of 27 common genes were obtained by further intersecting with IRGs and DEGs (Figure 4A).

Figure 4 Identification of potential diagnostic genes for sepsis. (A) Venn diagram of GSE69528, IMMPORT database and immune infiltrating cell related genes. (B) AUC value histogram of 10 key genes in 6 datasets. (C–H) ROC values of up-regulated genes (IL1R2, LTB4R, S100A11, S100A12, SORT1). (I–N) ROC values of down regulated genes (RASGRP1, CD3G, CD40LG, CD8A, PPP3CC).

Screening Key Genes According to the Diagnostic Value

To further screen the key genes with high diagnostic efficiency, we calculated the AUC values of 27 common genes in six datasets (GSE13904, GSE26440, GSE28750, GSE25504, GSE69528 and GSE145227). Finally, we identified ten genes with AUC values greater than 0.8 in six datasets and defined them as key genes (Figure 4B). Figure 4C–H shows the ROC curve of up-regulated genes (IL1R2, LTB4R, S100A11, S100A12 and SORT1); Figure 4I–N shows the ROC curve of down-regulated genes (RASGRP1, CD3G, CD40LG, CD8A and PPP3CC).

Subsequent Analysis of Key Genes: Expression, Regulation and Correlation with Immune Cells

In order to explore the expression in other data sets, upstream regulation mechanism and correlation with immune cells of these key genes, we made further analysis with R software. Among key genes, IL1R2, LTB4R, S100A11, S100A12 and SORT1 were up-regulated in sepsis, while RASGRP1, CD3G, CD40LG, CD8A and PPP3CC were down-regulated. Importantly, the differential expression patterns of these key genes were consistent in the above 6 datasets (Figure 5A). We studied TF-miRNA-mRNA interactions (Figure 5B) and detected central regulatory biomolecules (TFs and miRNAs) using topological parameters. Correlation analysis results show that there is a positive correlation between genes that are regulated in the same direction, and there is a significant negative correlation between up-regulated genes and down-regulated genes (Figure 5C). Correlation between key genes and immune infiltrating cells showed that the up-regulated genes were significantly positively correlated with macrophages, neutrophils and memory B cells, and negatively correlated with T cells (CD4/CD8) and NK cells, the results of down-regulated genes are just opposite to those of up-regulated genes (Figure 5D).

Figure 5 Expression and correlation analysis of 10 key genes. (A) Expression of 10 key genes in 6 datasets. (B) TF-miRNA-mRNA network construction (triangle: TF, square: miRNA, circle: mRNA). (C) Correlation analysis of 10 key genes. (D) Correlation between key genes and immune infiltrating cells, based on Pearson correlation analysis. Red nodes represent positive correlations, blue nodes represent negative correlations. *P < 0.05, **P < 0.01.

Functional Enrichment in Key Genes

We used AmiGO database to confirm GO term enrichment for identified key genes involved in sepsis (Table 2). To analyze the association between key genes and disease and GO term enrichment for key genes, we used CTD and amigo database for further research. The CTD databases show that key genes are mainly targeted in inflammatory and immune-related diseases, including sepsis, acute kidney injury, acute lung injury and multiple organ dysfunction (Table 3).

Table 2 Use the Amigo Database to Analyze the GO Terms of Keygenes

Table 3 Diseases Related to Key Gene Regulation

Predictive Value of Key Genes in Sepsis Survival

To determine the predictive value of key genes for sepsis survival, we performed univariate and multivariate logistic regression analysis. Among key genes, univariate regression analysis showed that 6 genes (LTB4R, S100A11, SORT1, RASGRP1, CD3G and CD8A) expression were significantly associated with sepsis survival. Multivariate regression analysis showed that S100A11 was an independent prognostic factor for sepsis survival (Table 4).

Table 4 Univariate and Multiple Regressive Analysis of Death for Sepsis

In the comparison of survival and non-survival groups of sepsis, the expression of down-regulated genes (RASGRP1, CD3G and CD8A) were significantly decreased in the survival group; the expression of up-regulated genes (LTB4R, S100A11, SORT1) increased significantly in the survival group (Figure 6).

Figure 6 (A-F) Expression of RASGRP1, CD8A, CD3G, SORT1, LTB4R, S100A11, respectively. *P < 0.05, **P < 0.01, ***P < 0.001.

Our results revealed that the 6 key genes expression levels are significantly associated with survival in sepsis. Most importantly, multivariate analysis showed that the S100A11 was an independent and favorable marker of prognosis in patients with sepsis. Therefore, the expression of S100A11 has the potential to predict the prognosis of patients with sepsis.

PCR Validation of Key Gene Expression Profiles

The expression of 10 key genes in sepsis patients (n = 6) and healthy controls (n = 6) was compared by real-time quantitative PCR. The results showed that compared with the control group, the expression of genes (RASGRP1, CD3G, CD40LG, CD8A and PPP3CC) was significantly down-regulated and genes (IL1R2, S100A12) was significantly up-regulated in patients with sepsis (P < 0.05; Figure 7). The expression results of these 7 key genes were consistent with the trend observed in microarray analysis.

Figure 7 (A-G) Expression of CD40LG, PPP3CC, RASGRP1, CD8A, CD3G, S100A12 and IL1R2, respectively. Notes: *P < 0.05, **P < 0.01, ***P < 0.001.

Discussion

Sepsis accounts for one-fifth of deaths worldwide, and the number of deaths from sepsis has doubled compared with estimates in previous years. Early diagnosis, prognostic evaluation and timely treatment of sepsis play an important role in reducing mortality. In essence, sepsis is a disorder of immune response caused by severe infection, but the specific mechanism is still not very clear. Despite years of research, the diagnosis and prediction of sepsis and septic shock still mainly rely on clinical criteria. Many biomarkers have been studied, none of which have the accuracy required for clinical practice.²⁵ In this study, bioinformatic analysis of sepsis-related gene expression data was performed to explore the function of immune cell infiltration in sepsis and to identify molecules involved in causing immune dysregulation in sepsis to assist in the early diagnosis of sepsis and to predict sepsis prognosis. We found that in patients with sepsis, there was an increased infiltration of intrinsic immune cells and a significant decrease in adaptive immune cell infiltration. We identified 10 key genes with high diagnostic efficacy associated with the immune infiltration characteristics of sepsis, and six of these may be key genes affecting the prognosis of sepsis.

By intersecting immune-related and differentially expressed genes, we identified immune-related differential genes associated with sepsis. Enrichment analysis revealed that these immune-associated differential genes are mainly involved in inflammatory and immunometabolic responses. In our analysis, immune cell infiltration in septic patients showed significantly higher levels of infiltration of certain intrinsic immune cells (eg, macrophages, neutrophils, mast cells) and reduced levels of infiltration of T cells (CD4, CD8) and NK cells. This is similar to the results of previous studies.²⁵ In the sepsis host, macrophages generate pro-inflammatory cytokines and chemokines which upregulate the inflammatory response and fight pathogens by generating a systemic inflammatory response in innate immune cells.²⁶ Neutrophils can migrate rapidly from the bone marrow in early sepsis, and within a few hours, the number of neutrophils can be 10 times higher than normal.²⁷ The number of circulating NK cells is significantly reduced in patients with sepsis, possibly due to increased apoptosis, and the results are associated with increased mortality.²⁸ It has been found that reduced and dysfunctional NK cells are closely associated with sepsis-induced immunosuppression, leading to a greater susceptibility to infection and a poorer prognosis for patients.^29,30 It is evident that these intrinsic immune cells play an important role in the regulation of inflammation in sepsis. On the other hand, sepsis increases apoptosis, which may lead to a decrease in the number of CD4 T cells,³¹ and significantly reduces the proliferative capacity of CD4 T cells.³² Also, sepsis-induced apoptosis reduces the number of naive and memory CD8 T cells.^33,34 Condotta et al observed increased expression of co-suppressor molecules (eg, PD-1, LAG-3 and 2B4) and reduced release of effector cytokines (eg, IFN-γ, TNF-α and IL-2), suggesting that CD8 T cells are depleted in mice and sepsis patients.^35–37 Overall, our findings are largely consistent with previous reports that sepsis leads to a significant reduction in CD4 and CD8 T cells, which is a major driver of impaired immune responses and reduced host antimicrobial capacity.^38–40

We further identified 10 key genes (IL1R2, LTB4R, S100A11, S100A12, SORT1, RASGRP1, CD3G, CD40LG, CD8A and PPP3CC) with high diagnostic efficacy by AUC analysis, all of which had AUC values greater than 0.8 for the diagnosis of sepsis in six datasets, and therefore they may be early diagnostic markers for sepsis. Analysis of the correlation between these 10 key genes and immune infiltrating cells revealed that their differential expression was mainly associated with high levels of infiltration of macrophages and neutrophils and low levels of infiltration of T cells (CD4, CD8) and NK cells significantly. We speculate that upregulation of the key genes IL1R2, LTB4R, S100A11, S100A12, SORT1 and downregulation of the key genes RASGRP1, CD3G, CD40LG, CD8A and PPP3CC may play an important role in the enhancement of intrinsic immunity and suppression of adaptive immunity during the development of sepsis. Specifically, our PCR results showed that the expression levels of 7 key genes (IL1R2, S100A12, RASGRP1, CD3G, CD40LG, CD8A and PPP3CC) were consistent with the results of microarray analysis, while other three genes (LTB4R, S100A11, SORT1) were not significantly different. This inconsistency may be related to different assay methods, sample size, patient heterogeneity and disease duration. Later experiments with larger sample sizes are needed for validation.

The high mortality rate of sepsis has long plagued the world. Immunosuppression is now largely recognised as a major cause of sepsis-related death, with T-cell depletion being one of the most serious factors.^5,41 Significant reductions in NK-cell numbers have also been associated with increased mortality.²⁸ However, the exact mechanisms are not well understood. To further understand the association between key genes and the survival prognosis of sepsis, we used logistic regression analysis to identify 10 key genes, and the results suggested that six key genes (LTB4R, S100A11, SORT1, RASGRP1, CD3G and CD8A) might affect the survival prognosis of sepsis, among which S100A11 might be an independent prognostic factor for survival. These six key genes were significantly differentially expressed in the survival group compared to the non-survival group of sepsis, further suggesting that they may be prognostic markers of sepsis.

SA10012 is mainly expressed and secreted by neutrophils and is an important pro-inflammatory factor that plays an important role in microbial resistance, regulation of apoptosis and involvement in inflammatory and immune responses.⁴² It has been shown that S100A12 expression levels are significantly upregulated in neonatal septic rats and that it has a positive pro-inflammatory role in the inflammatory response, a role that is important in the development and progression of sepsis.⁴² Whereas S100A11 was mainly involved in the process of tumor development in previous studies, it is not clear whether it is associated with disease regulation in sepsis. We found that S100A11 expression was upregulated in septic patients and independently affected the prognosis of septic patients, the mechanism of which needs further investigation. RASGRP1 is a key signaling molecule essential for lymphocyte function. Immunologically, RASGRP1 deficiency is characterized by reduced CD4 and CD8 T-cell lymphocytes and impaired proliferation, as well as impaired B-cell development and proliferation, and abnormal NK cytotoxic function.^43,44 RASGRP1-deficient mice have shown significant T-cell lymphopenia and immunodeficiency⁴⁵ and are more susceptible than wild-type mice to sepsis infections caused by cecum ligation and puncture.⁴⁶ Similarly, we found reduced levels of RASGRP1 expression in sepsis patients, suggesting that impaired immunity in sepsis patients may be associated with low RASGRP1 expression. Therefore, we speculate that RASGRP1 may play a very critical role in the development of sepsis. LTB4R is the receptor for LTB 4, which induces the recruitment of CD4 and CD8 effector T cells as well as helper T cells to inflammatory tissues,^47,48 and the recruited immune cells will promote inflammation by secreting cytokines and chemokines. Sort1 may participate in innate immunity by regulating inflammation and phagocytosis, as well as participating in adaptive immunity by controlling the maturation of immune cells and regulating the activation of T cells and NK cells.⁴⁹ CD3G is involved in the formation of T cell receptor CD3 complex and is an important regulator of T cell development and differentiation.⁵⁰ Its defects are related to T cell immunodeficiency.⁵¹ CD8A is primarily involved in cell-mediated immune defence and T cell development.^52,53 The role of these key genes in sepsis is rarely studied, but the results of our analysis suggest that they may contribute to the early diagnosis and prognostic assessment of sepsis.

Our research also has some limitations. First of all, because our research is based on multiple public datasets, the course of sepsis may be inconsistent. Secondly, due to the small sample size and short study time, there is a lack of patient data, such as blood culture results and mortality. We will increase sample size, extend research time, and add more patient data in subsequent research.

Conclusions

In conclusion, we identified ten key genes (IL1R2, LTB4R, S100A11, S100A12, SORT1, RASGRP1, CD3G, CD40LG, CD8A and PPP3CC) as potential diagnostic biomarkers of sepsis. Six key genes (LTB4R, S100A11, SORT1, RASGRP1, CD3G and CD8A) might affect the survival prognosis of sepsis, among which S100A11 might be an independent prognostic biomarker for survival. These genes may aid in the early diagnosis and assessment of sepsis, and could serve as interesting targets for further in-depth research.