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

Purpose: Sepsis is a serious life-threatening condition characterised by multi-organ failure due to a disturbed immune response caused by severe infection. The pathogenesis of sepsis is unclear. The aim of this article is to identify potential diagnostic and prognostic biomarkers of sepsis to improve the survival of patients with sepsis.
Methods: We downloaded 7 datasets from Gene Expression Omnibus database and screened the immune-related differential genes (IRDEGs). The related functions of IRDEGs were analyzed through Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). CIBERSORT was used to evaluate the infiltration of the immune cells, and Pearson algorithm of R software was used to calculate the correlation between the immune cell content and gene expression to screen the genes most related to immune cells in sepsis group, which were intersected with IRDEGs to obtain common genes. Key genes were identified from common genes based on the area under the receiver operating characteristic curve (AUC) greater than 0.8 in the 6 datasets. We then analyzed the predictive value of key genes in sepsis survival. Finally, we verified the expression of key genes in patients with sepsis by PCR analysis.
Results: A total of 164 IRDEGs were obtained, which were associated mainly with inflammatory and immunometabolic responses. Ten key genes (IL1R2, LTB4R, S100A11, S100A12, SORT1, RASGRP1, CD3G, CD40LG, CD8A and PPP3CC) were identified with high diagnostic efficacy. Logistic regression analysis revealed that six of the key genes (LTB4R, S100A11, SORT1, RASGRP1, CD3G and CD8A) may affect the survival prognosis of sepsis. PCR analysis confirmed that the expression of seven key genes (IL1R2, S100A12, RASGRP1, CD3G, CD40LG, CD8A and PPP3CC) was consistent with microarray outcome.
Conclusion: This study explored the immune and metabolic response mechanisms associated with sepsis, and identified ten potential diagnostic and six prognostic genes.

Introduction

Sepsis is a common clinical emergency and critical illness, and is regarded as one of the life-threatening infectious diseases.1 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.2

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.3 In other terms, sepsis progresses from excessive immune activation to widespread immunosuppression.4 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.3 Growing evidence supports immunosuppression as a major cause of sepsis-related deaths,5 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.9 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.12 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” package13 and volcano plots of DEGs were plotted using the “ggplot2” package14 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 (), 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.15 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 R16 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 TarBase17 and mirtarbase18 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.19 All of the above use NetworkAnalyst to analyze topology parameters (ie, degree and medium centrality).20 We further used Cytoscape software21 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; ) is a public database describing gene–disease interactions and providing information on gene–disease interactions.22 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; ) is a web application that is primarily used to query, browse, visualize and download gene ontologies and annotations.23 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,24 (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

PCR Primer	Primer Sequence	Product Length (bp)
h-IL1R2-F	GCGGGAGTTCAGGCTGGAAG	184
h-IL1R2-R	TGGCAGAAGCCACAGAGCAC
h-LTB4R-F	GCCCAAGGCACCTGGAGTTT	162
h-LTB4R-R	CGCCTTGGTGCGTAGCTTCT
h-S100A11-F	GACTGAGCGGTGCATCGAGT	162
h-S100A11-R	ATGCGGTCAAGGACACCAGG
h-S100A12-F	CGGAAGGGGCATTTTGACACC	187
h-S100A12-R	TCAGCGCAATGGCTACCAGG
h-SORT1-F	GGCTCCTGCAAAGCTGACCT	200
h-SORT1-R	TGGGCCATGCTCCATGTGTC
h-RASGRP1-F	GGCCTGGGCTTTCCTCACAA	188
h-RASGRP1-R	ACTGGGTTCTTGGCTCGCTT
h-CD3G-F	GGAGTTCGCCAGTCGAGAGC	185
h-CD3G-R	AGGAGGAGAACACCTGGACTAC
h-CD40LG-F	GCAACCAGGTGCTTCGGTGT	168
h-CD40LG-R	GTGGGCTTAACCGCTGTGCT
h-CD8A-F	GGAGGCACCCATGCCATCTC	190
h-CD8A-R	TGAGACAGGGGCCTCGGAAA
h-PPP3CC-F	ACACTGTCCGAGGGTGCTCT	144
h-PPP3CC-R	AGCCTGTGGCTTGGCTCTTC
h-GAPDH-F	GACCTGACCTGCCGTCTA	83
h-GAPDH-R	AGGAGTGGGTGTCGCTGT

Abbreviation: PCR, polymerase chain reaction.

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

Gene	GO Class (Direct)	Evidence	Reference
IL1R2	Interleukin-1 receptor activity	IBA	PMID:21873635
LTB4R2	Neuropeptide signaling pathway	IBA	PMID:21873635
	Leukotriene B4 receptor activity	IBA	PMID:21873635
	G protein-coupled peptide receptor activity	IBA	PMID:21873635
	Integral component of plasma membrane	IBA	PMID:21873635
	Inflammatory response	IBA	PMID:21873635
S100A11	Calcium-dependent protein binding	IBA	PMID:21873635
	Extracellular space	IBA	PMID:21873635
	S100 protein binding	IBA	PMID:21873635
	Calcium ion binding	IBA	PMID:21873635
	Cytoplasm	IBA	PMID:21873635
	Calcium ion binding	IBA	ZFIN:ZDB-PUB-020724-1
	Transition metal ion binding	IBA	ZFIN:ZDB-PUB-020724-1
	Extracellular space	IBA	PMID:21873635
	Calcium-dependent protein binding	IBA	PMID:21873635
	Perinuclear region of cytoplasm	IBA	PMID:21873635
S100A12	Calcium ion binding	IEA	GO_REF:0000002
	RAGE receptor binding	IEA	GO_REF:0000107
	Xenobiotic metabolic process	IEA	GO_REF:0000107
	Killing of cells of other organism	IEA	GO_REF:0000107
	Positive regulation of I-kappaB kinase/NF-kappaB signaling	IEA	GO_REF:0000107
	Defense response to fungus	IEA	GO_REF:0000107
	Antimicrobial humoral immune response mediated by antimicrobial peptide	IEA	GO_REF:0000107
	Nucleus	IEA	GO_REF:0000107
	Cytoplasm	IEA	GO_REF:0000107
	Protein binding	IPI	PMID:10399917
SORT1	Golgi apparatus	IBA	PMID:21873635
	Post-Golgi vesicle-mediated transport	IBA	PMID:21873635
	Integral component of membrane	IBA	PMID:21873635
	Vesicle organization	IBA	PMID:21873635
	Endocytosis	IBA	PMID:21873635
	Golgi to endosome transport	IBA	PMID:21873635
	Golgi to endosome transport	IBA	PMID:21873635
	Integral component of membrane	IBA	PMID:21873635
	Post-Golgi vesicle-mediated transport	IBA	PMID:21873635
	Vesicle organization	IBA	PMID:21873635
RASGRP1	Guanyl-nucleotide exchange factor activity	IMP	PMID:23908768
	Calcium ion binding	IMP	PMID:23908768
	Zinc ion binding	IMP	PMID:23908768
	Lipid binding	TAS	PMID:9582122
	Diacylglycerol binding	IMP	PMID:23908768
	Phosphatidylcholine binding	IMP	PMID:23908768
	Identical protein binding	IPI	PMID:23908768
CD3G	Alpha-beta T cell receptor complex	IBA	PMID:21873635
	Positive thymic T cell selection	IBA	PMID:21873635
	Cell surface receptor signaling pathway	IBA	PMID:21873635
	External side of plasma membrane	IBA	PMID:21873635
	Transmembrane signaling receptor activity	IBA	PMID:21873635
	Positive thymic T cell selection	IBA	PMID:21873635
	Alpha-beta T cell receptor complex	IBA	PMID:21873635
	Transmembrane signaling receptor activity	IBA	PMID:21873635
	External side of plasma membrane	IBA	PMID:21873635
	Cell surface receptor signaling pathway	IBA	PMID:21873635
CD40LG	Membrane	IEA	ZFIN:ZDB-PUB-020724-1
	Cytokine activity	IEA	ZFIN:ZDB-PUB-020723-1
	Tumor necrosis factor receptor binding	IEA	ZFIN:ZDB-PUB-020724-1
	Immune response	IEA	ZFIN:ZDB-PUB-020724-1
	Integral component of membrane	IEA	ZFIN:ZDB-PUB-020723-1
	Membrane	IEA	ZFIN:ZDB-PUB-020723-1
	Extracellular space	IEA	ZFIN:ZDB-PUB-020723-1
	Positive regulation of immunoglobulin production	IDA	PMID:19494299
	CD40 receptor binding	IDA	PMID:19494299
	External side of plasma membrane	IDA	PMID:19494299
CD8A	Plasma membrane	IEA	ZFIN:ZDB-PUB-120306-4
	Integral component of membrane	IEA	ZFIN:ZDB-PUB-020723-1
	Membrane	IEA	ZFIN:ZDB-PUB-020723-1
	Molecular function	ND	ZFIN:ZDB-PUB-031118-1
	MHC class I protein complex binding	IEA	GO_REF:0000107
	External side of plasma membrane	IEA	GO_REF:0000107
	Integral component of membrane	IEA	GO_REF:0000043
	Plasma membrane raft	IEA	GO_REF:0000107
	Integral component of membrane	IEA	GO_REF:0000043
	MHC class I protein complex binding	IEA	GO_REF:0000107
PPP3CC	Cytoplasm	IBA	PMID:21873635
	Calcineurin-mediated signaling	IBA	PMID:21873635
	Calmodulin binding	IBA	PMID:21873635
	Calmodulin-dependent protein phosphatase activity	IBA	PMID:21873635
	Calcineurin complex	IBA	PMID:21873635
	Calmodulin binding	IEA	ZFIN:ZDB-PUB-020723-1
	Phosphoprotein phosphatase activity	IEA	ZFIN:ZDB-PUB-020723-1
	Hydrolase activity	IEA	ZFIN:ZDB-PUB-020723-1
	Calmodulin-dependent protein phosphatase activity	IEA	ZFIN:ZDB-PUB-020723-1
	Hydrolase activity	IEA	ZFIN:ZDB-PUB-020723-1

Abbreviation: GO, Gene Ontology.

Table 3

Diseases Related to Key Gene Regulation

Gene	Disease Name	Disease ID	Inference Score	Reference Count
IL1R2	Inflammation	MESH:D007249	143.71	259
	Immune System Diseases	MESH:D007154	57.41	23
	Sepsis	MESH:D018805	14.49	18
	Immunologic Deficiency Syndromes	MESH:D007153	36.92	6
	Autoimmune Diseases	MESH:D001327	32.27	12
RASGRP1	Acute Kidney Injury	MESH:D058186	97.07	391
	Acute Lung Injury	MESH:D055371	52.73	68
	Sepsis	MESH:D018805	27.87	30
	Immune System Diseases	MESH:D007154	41.24	17
	Anti-Neutrophil Cytoplasmic Antibody-Associated Vasculitis	MESH:D056648	12.88	51
CD3G	Inflammation	MESH:D007249	127.93	209
	Acute Kidney Injury	MESH:D058186	63.38	183
	Lung Injury	MESH:D055370	74.6	77
	Sepsis	MESH:D018805	15.26	17
	Autoimmune Diseases	MESH:D001327	26.01	8
CD8A	Acute Kidney Injury	MESH:D058186	60.12	269
	Inflammation	MESH:D007249	122.47	265
	Acute Lung Injury	MESH:D055371	55.32	39
	Sepsis	MESH:D018805	33.58	21
	Immune System Diseases	MESH:D007154	37.71	15
CD40LG	Inflammation	MESH:D007249	176.78	321
	Acute Kidney Injury	MESH:D058186	101.01	193
	Acute Lung Injury	MESH:D055371	55.62	70
	Sepsis	MESH:D018805	43.35	37
	Immune System Diseases	MESH:D007154	40.23	13
LTB4R	Inflammation	MESH:D007249	66.35	105
	Acute Kidney Injury	MESH:D058186	27.17	120
	Liver Failure, Acute	MESH:D017114	19.26	160
	Acute Lung Injury	MESH:D055371	16.25	8
	Sepsis	MESH:D018805	2.6	2
PPP3CC	Acute Kidney Injury	MESH:D058186	44.08	235
	Liver Failure, Acute	MESH:D017114	29.48	167
	Inflammation	MESH:D007249	91.18	136
	Immune System Diseases	MESH:D007154	29.59	15
	Sepsis	MESH:D018805	15.16	9
S100A11	Acute Kidney Injury	MESH:D058186	150.56	420
	Inflammation	MESH:D007249	237.55	337
	Acute Lung Injury	MESH:D055371	55.04	101
	Sepsis	MESH:D018805	39.53	38
	Immune System Diseases	MESH:D007154	46.9	20
S100A12	Inflammation	MESH:D007249	58.08	173
	Acute Kidney Injury	MESH:D058186	20.05	106
	Acute Lung Injury	MESH:D055371	20.15	44
	Shock, Septic	MESH:D012772	3.32	15
	Sepsis	MESH:D018805	6.66	13
SORT1	Inflammation	MESH:D007249	173.13	229
	Acute Kidney Injury	MESH:D058186	91.31	214
	Acute Lung Injury	MESH:D055371	46.1	70
	Sepsis	MESH:D018805	26.53	17
	Multiple Organ Failure	MESH:D009102	33.14	8

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

Characteristics	Total (N)	Univariate Analysis		Multivariate Analysis
		Hazard Ratio (95% CI)	P value	Hazard Ratio (95% CI)	P value
Sex	127
Female	75	Reference
Male	52	2.376 (1.092–5.167)	0.029	3.766 (1.473–9.625)	0.006
APACHE II score	127	1.121 (1.047–1.201)	0.001	1.098 (1.023–1.179)	0.009
Neutrophil proportion	126	1.350 (0.055–33.301)	0.854
IL1R2	127	0.704 (0.208–2.389)	0.574
LTB4R	127	0.376 (0.187–0.759)	0.006	1.639 (0.460–5.830)	0.446
S100A11	127	0.384 (0.217–0.678)	<0.001	0.329 (0.120–0.898)	0.030
S100A12	127	0.723 (0.473–1.105)	0.134
SORT1	127	0.408 (0.223–0.749)	0.004	0.668 (0.300–1.489)	0.324
RASGRP1	127	3.361 (1.107–10.199)	0.032	1.114 (0.117–10.569)	0.925
CD3G	127	6.421 (2.150–19.180)	<0.001	2.487 (0.521–11.876)	0.254
CD40LG	127	0.410 (0.035–4.817)	0.478
CD8A	127	1.606 (1.076–2.397)	0.020	0.982 (0.444–2.175)	0.965
PPP3CC	127	2.186 (0.476–10.035)	0.315

Notes: APACHE II score: Acute Physiology and Chronic Health Evaluation II score. All P values less than 0.05 are shown in bold font.

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.

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.25 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.25 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.26 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.27 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.28 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,31 and significantly reduces the proliferative capacity of CD4 T cells.32 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.28 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.42 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.42 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 immunodeficiency45 and are more susceptible than wild-type mice to sepsis infections caused by cecum ligation and puncture.46 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.49 CD3G is involved in the formation of T cell receptor CD3 complex and is an important regulator of T cell development and differentiation.50 Its defects are related to T cell immunodeficiency.51 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.