Screening of Key Genes of Sepsis and Septic Shock Using Bioinformatics Analysis

Introduction

Sepsis is defined as life-threatening organ dysfunction caused by a dysregulated host response to infection.¹ In 2017, approximately 48.9 million sepsis cases and 11.0 million sepsis-related deaths were recorded worldwide, accounting for 19.7% of all global deaths.² According to a recently published systematic review, sepsis with organ dysfunction occurring during hospitalization affects 24.4% of all patients in intensive care units and neonatal intensive care units.³ Septic shock is defined as sepsis complicated by hyperlactatemia and simultaneous hypotension requiring vasopressin therapy. The in-hospital mortality rate of septic shock is 30–50%.⁴

The World Health Organization states that the annual number of deaths due to sepsis worldwide is approximately 6 million, most of which are preventable.⁵ The identification of altered molecular signatures and biochemical pathways in sepsis patients has driven interest in the discovery of novel biomarkers.⁶ Early recognition of sepsis is essential for timely treatment and can improve sepsis outcomes.^7,8 Mortality risk is higher when sepsis treatment is delayed.⁹ Analyzing the molecular characteristics of patients is an effective method to screen for potential diagnostic and prognostic targets.¹⁰

Patients with sepsis often present acquired immunodeficiency.¹¹ Significant lymphopenia has been reported in sepsis patients, and apoptosis may be a major mechanism of lymphocyte death in sepsis.¹² A number of anti-inflammatory and immunostimulatory agents are currently in clinical trials for sepsis and septic shock.^13,14 At the same time, sepsis may result in prolonged methylation changes in monocytes, which may help stabilize them in an immunotolerant state.¹⁵ Multiple studies have shown that epigenetic alterations in sepsis may affect prognosis and serve as diagnostic biomarkers.^16,17

Improved infection prevention and control strategies are urgently needed to reduce the health care burden associated with the development and progression of sepsis and septic shock, as are approaches to early diagnosis and appropriate treatment. The present study used bioinformatics to explore biomarkers and potential therapeutic targets associated with the development of sepsis.

Materials and Methods

Sepsis Data Collection

Data were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The GSE13904 dataset included array-based gene expression profiles of whole blood from 18 normal children, 52 children with sepsis, and 106 children with septic shock. The GSE25504 dataset included array-based gene expression profiles of whole blood from 37 control and 25 infected human neonates. The GSE9960 dataset included array-based gene expression profiles of peripheral blood mononuclear cells from 54 adult sepsis patients and 16 controls. The raw data in these three datasets were processed using the lumi package in R.¹⁸ The GSE154918 dataset included gene expression profiles based on high-throughput sequencing of peripheral blood from 65 adult sepsis patients and 40 controls. Raw counts were pre-processed and normalized using DESeq2.¹⁹

Coexpression Analysis

The variance of gene expression in GSE13904 (children) was first calculated, and the top 20% of the greatest variance was screened to construct a coexpression network through weighted gene coexpression network analysis (WGCNA) using the WGCNA package in R.²⁰ The connections among different pairs of genes were identified and weighted based on the correlated expression levels across patients and control samples. Then the adjacency matrix was converted into a topological overlap matrix (TOM) to detect gene connectivity in the network. Genes were separated into different clusters (modules) based on their connectivity and covariance coefficients, then hierarchically clustered. The WGCNA package was employed to test the independence and average connectivity of different modules under different power values, and the power values corresponding to an independence index of R² = 0.9 were selected. The minimum size of the gene dendrogram was 30. Potential correlations between modules and clinical characteristics of patients were explored through Pearson correlation analysis.

Difference Analysis

Gene expression profiles were constructed for each sample of septic shock patients, sepsis patients and controls. Differential analysis between sepsis patients and healthy individuals in GSE13904 (children) or GSE154918 (adults) was performed after using the limma package in R.²¹ A screening threshold of P < 0.05 was set to obtain the differentially expressed genes (DEGs) between sepsis patients and controls. Then, DEGs whose expression differed in the same direction in patients (up- or down-regulation) were screened for intersection analysis with module genes to obtain the common genes. DEGs between septic shock patients and sepsis patients in GSE13904 (children) was obtained using the limma package in R and a significance criterion of P < 0.05.

Enrichment Analysis

We used the clusterProfiler package in R²² to perform gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of common genes. P < 0.05 was considered statistically significant for enrichment. The significantly enriched terms were subjected to gene set variation analysis (GSVA)²³ to calculate the score for the enrichment using the gene expression profile. The up-regulation or down-regulation of enrichment terms in septic patients relative to controls was evaluated using the GSVA score in the limma package in R. Gene set enrichment analysis (GSEA) was used to identify KEGG pathway using GSEA software. The fgsea package in R was used to display the results of the GSEA.

Single Sample Gene Set Enrichment Analysis (ssGSEA)

The marker gene sets for immune cell types were obtained from Bindea et al²⁴. The gene signatures expressed by the immune cell populations were quantified using the single-sample gene set enrichment analysis (ssGSEA) function in GSVA. Differences in immune cell infiltration between sepsis patients and healthy controls were then calculated, and a significance criterion of P < 0.05 was applied. Correlations between levels of infiltration by different types of immune cells were calculated using Pearson’s correlation, based on ssGSEA scores for the immune cell types.

PPI Network Analysis

The common genes were analyzed using the online tool STRING (https://string-db.org), and the protein-protein interaction (PPI) network based on the criterion of a combined score > 900 was constructed and displayed using Cytoscape software. Genes were ranked based on their degree of connectivity with other genes. The pROC package in R²⁵ was used to calculate the the area under the receiver operating characteristic curve (AUC) for the top 100 genes showing the greatest connectivity in the PPI network.

Short Time-Series Expression Miner (STEM) Analysis

We performed STEM analysis to cluster the common genes in healthy individuals, sepsis patients and septic shock patients in the GSE13904 database (children).²⁶ Significant clustering was defined as P < 0.05. The significantly clustered genes showed a trend of gradual up- or down-regulation as one moved from controls to sepsis patients and then to septic shock patients.

Methylation Analysis

The GSE138074 dataset contains array-based gene methylation profiles of monocytes from 14 adult septic patients and 11 healthy donors. Data were processed using the ChAMP package in R.²⁷ Differences in methylation levels between sepsis patients and healthy donors were calculated using the limma package in R.

Results

Coexpression Network of Sepsis-Related Genes

The flow chart of this study is shown in Figure 1. We calculated variance of gene expression in GSE13904 (children) and then screened 4093 genes in the top 20% greatest variance genes in order to identify coexpression patterns. WGCNA identified a coexpression network containing 2550 genes (Figure 2A). These genes clustered into eight coexpression modules: the brown module had the strongest positive correlation with sepsis shock, whereas the pink module had the strongest negative correlation with sepsis shock (Figure 2B). By calculating the expression profiles of the modules in different clinical samples, we identified the expression trend for each module in sepsis development (Figure 2C). The brown module showed gradual up-regulation, while the pink module showed gradual down-regulation, in the order: controls < patients with sepsis < patients with septic shock.

Figure 1 Flowchart of this study. The following datasets were used for the identification of potential diagnostic genes and mechanisms associated with the development of sepsis: GSE13904 (children), GSE25504 (children), GSE9960 (adults) and GSE154918 (adults). Abbreviations: AUC, area under the receiver operating characteristic curve; GSVA, gene set variation analysis; PPI, protein-protein interaction; STEM, short time-series expression miner; WGCNA, weighted gene coexpression network analysis.

Figure 2 WGCNA of genes showing large expression variance in sepsis patients. (A) Network heatmap of unsupervised cluster analysis of module genes. Different colors in columns and rows represent different modules. (B) Correlation of modules with clinical traits. Each column represents a different module; each row, a different clinical phenotype. Red indicates positive correlation; blue, negative correlation. (C) Progressive up- or down-regulation in a module in the trend: healthy controls < sepsis < sepsis shock.

DEGs in Sepsis

Of the four sets of sepsis data that we obtained, GSE13904 (children) and GSE154918 (adults) had the largest samples (Figure 3A). In GSE13904 (children), we obtained 7850 DEGs between sepsis and control samples (Figure 3B, Table S1). In GSE154918 (adults), we obtained 12496 DEGs between sepsis and control samples (Figure 3B, Table S2). By comparing the DEGs between adult and pediatric sepsis, we found 5638 genes that were common to both, 2212 that may be specific to the pediatric condition and 6857 that may be specific to the adult condition (Figure 3C). Genes at the intersection may be strongly associated with sepsis in adults and children. In addition, we screened 5143 DEGs that were simultaneously up- or down-regulated in GSE13904 (children) and GSE154918 (adults). Intersection analysis of these genes with module genes yielded 1274 common genes for subsequent analysis (Figure 3D).

Figure 3 Identification of common genes in sepsis. (A) Petal plot of sample size for sepsis and control groups in four datasets. Red petals represent healthy controls; blue petals, sepsis patients. (B) Manhattan plot of differentially expressed genes between sepsis patients and healthy controls in GSE13904 and GSE154918. Up, up-regulation in sepsis; Down, down-regulation in sepsis; Not, no significant difference from controls. The three genes in each group showing the largest fold difference are labeled. (C) Venn diagram of differentially expressed genes in GSE13904 (children) and GSE154918 (adults). (D) Venn diagram showing DEGs in GSE13904 (children) and GSE154918 (adults) as well as module genes, that are common to adults and children or specific to each one of those groups.

Biological Functions of Selected Genes

Enrichment analysis showed that, based on the common genes identified above, defensive responses to bacteria, complement-dependent cytotoxicity, and canonical glycolysis were significantly up-regulated in sepsis patients, whereas interleukin-17 production, B cell activation, and T cell receptor signaling were significantly down-regulated (Figure 4A). KEGG results showed that Hypoxia-Inducible Factor 1 (HIF-1) signaling pathway, Tumor Necrosis Factor (TNF) signaling pathway, and glycolysis/gluconeogenesis were significantly up-regulated in sepsis, whereas Th1 and Th2 cell differentiation, Th17 cell differentiation, and T cell receptor signaling pathways were significantly down-regulated (Figure 4B). KEGG results in GSEA similarly found that metabolism-related pathways were up-regulated while immune related pathways were down-regulated (Figure 4C).

Figure 4 Biological functions and KEGG pathways enriched in common genes. (A) Common genes of sepsis patients relative to controls involved in up- or down-regulated biological processes, as quantified by gene set variation analysis (GSVA). FC, fold change. (B) Common genes of sepsis patients relative to controls involved in up- or down-regulated KEGG pathways, as quantified by single-sample GSVA. FC, fold change. (C) Common genes involved in up- or down-regulated KEGG pathways of GSEA results in sepsis patients relative to controls. P < 0.05 was considered statistically significant.

Immune Cell Infiltration in Sepsis

In the enrichment results, we found that a large number of immune-related biological functions were significantly down-regulated in sepsis samples. Therefore, we assessed immune cell infiltration in sepsis patients (Figure 5A). Unsurprisingly, sepsis was associated with significantly reduced infiltration by most types of immune cells, except that infiltration by macrophages, mast cells, and neutrophils was elevated. Cytotoxic T cell types were defined based on gene markers, and CD8+ T cells were treated as one subtype of T cells with different gene markers than other subtypes.²⁴ The ssGSEA scores for these different subtypes were used to generate immune cell interaction networks. We grouped the differentially infiltrated immune cells into four clusters (Figure 5B).

Figure 5 Immune cell suppression in sepsis patients. (A) Differences in immune cell infiltration between sepsis and controls. Each row represents a type of immune cell; each column, a different dataset. Red nodes represent significant up-regulation in sepsis; blue nodes, significant down-regulation. (B) Clustering and correlation of immune cell types based on infiltration levels. The size of each node represents the single-sample gene set enrichment analysis (ssGSEA) scores of each immune cell type, transformed by log10 (Log rank test P-value). Connections between immune cell types represent interactions between the two. The thickness of the line indicates the strength of the correlation, based on Pearson correlation analysis. Red lines represent positive correlations; blue lines, negative correlations. Immune cell cluster A is shown in yellow; cell cluster B, blue; cell cluster C, red; and cell cluster D, brown.

PPI Network of Common Genes

In addition, we performed PPI network analysis of the common genes and ranked the connectivity degree between the genes in the network (Figure 6A). We screened the top 100 genes with the greatest connectivity (Figure 6B). By calculating the AUC values of these genes in GSE13904 (children) and GSE154918 (adults), we identified genes with AUC values greater than 0.9 in both sets of data, which we defined as key genes (Figure 6C). Among them, MAPK14, FGR, RHOG were up-regulated in sepsis, while LAT, PRKACB, UBE2Q2, ITK, IL2RB, and CD247 were down-regulated (Figure 6D). Importantly, the patterns of differential expression of these key genes were consistent across GSE13904, GSE154918, GSE25504, and GSE9960. Correlation analysis found that up-regulated genes had stronger negative correlation with immune cells, while down-regulated genes had stronger positive correlation with immune cells, especially T cells (Figure 6E).

Figure 6 Identification of potential diagnostic genes for sepsis. (A) PPI network of common genes (top 100 with highest degree of connection). Darker color indicates higher degree of connectivity. (B) Heatmap of the top 100 genes based on degree of connectivity in the PPI network. Data from sepsis patients are shown in yellow; sepsis shock patients, green; and controls, brown. (C) Genes with areas under the receiver operating characteristic curve (AUCs) greater than 0.9 in GSE13904 (children) and GSE154918 (adults). Red represents genes with up-regulated expression in sepsis; blue, genes with down-regulated expression. The lengths of the bars represent the mean AUC values of the genes. (D) Differential expression of key genes between sepsis patients and controls. Red nodes represent genes with up-regulated expression in sepsis; blue nodes, genes with down-regulated expression; grey nodes, no significant difference from controls. (E) Correlation between key genes and immune infiltrating cells, based on Pearson correlation analysis. Red nodes represent positive correlations; blue nodes, negative correlations. *P < 0.05, **P < 0.01.

Genes Associated with Sepsis Development

To identify the change in gene expression during the course of development from sepsis to septic shock, we utilized STEM software to identify 712 common genes showing progressive dysregulation during sepsis progression. These genes fell into five significant clusters (Figure 7A). Interestingly, all the key genes that we identified above were detected in these clusters (Figure 7B). The expression of MAPK14, FGR, RHOG was gradually up-regulated from control to sepsis then to septic shock, while the expression of LAT, PRKACB, UBE2Q2, ITK, IL2RB, and CD247 was gradually down-regulated (Figure 7C). The pathway heatmap built by GSEA algorithm showed that starch and sucrose metabolism, complement and coagulation cascades, and legionellosis were gradually up-regulated during the development of sepsis, while Th1 and Th2 cell differentiation, IgA production in the intestinal immune network, as well as antigen processing and presentation were gradually down-regulated (Figure 7D).

Figure 7 Identification of persistently dysregulated genes during progression from sepsis to septic shock. (A) Heatmap of persistently dysregulated genes identified by STEM in the trend: healthy controls < sepsis patients < septic shock patients. Gene sets were arranged based on cluster assignment in order to generate simplified expression profiles. We graphically depicted only five clusters with >10 genes. (B) The box plots of STEM genes in five clusters. Line plots and box plots were used to display log2 (fold change). Representative genes were highlighted using red lines. Key genes were in significant STEM clustering and located on the right side of the box map. *P < 0.05, **P < 0.01. (C) Expression of key genes were persistently elevated during sepsis development. (D) Signaling pathways persistently up- or down-regulated as sepsis develops. Red in the heatmap represents an enrichment score greater than 0; blue represents an enrichment score less than 0.

Methylation Marks in Sepsis

To identify the sepsis key genes whose methylation may be modified in the disease, we first analyzed differentially methylated positions (DMPs) between sepsis samples and controls in GSE138074 (adult) (Figure 8A). DMPs with delta beta values in the opposite direction to the expression value of common DEGs were filtered and identified as methylation marks. A total of 1313 methylation marks were identified (Figure 8B). These included MAPK14, FGR, and CD247. The methylation levels of MAPK14 (cg18213931) and FGR (cg16922167) were lower in sepsis than control, and that of CD247 (cg21161394) was higher in sepsis (Figure 8C).

Figure 8 Screening of key methylation markers for sepsis in GSE138074. (A) Differentially methylated positions (DMPs) between sepsis and controls. Of all DMPs, 49.58% were hypermethylated and 50.42% were hypomethylated. The 10 sites with the highest methylation levels are labeled. Each feature color represents a different methylation site. (B) The methylation levels and the expression levels of methylation marks. Key genes are labeled. (C) Differences in key gene methylation levels between sepsis and controls. P for CD247 was 0.039; MAPK14, 0.047; FGR, 0.047.

Discussion

Early recognition of sepsis is key to timely treatment. In this study, we subjected sepsis-related sequencing data to bioinformatics analysis in order to explore the sepsis-related molecular dysregulation mechanisms. Compared to other studies using the same datasets,^28,29 the present work used WGCNA and PPI networks to identify potential diagnostic markers for sepsis not only in adults but also in children. We further identified potential biomarkers associated with septic shock as well as methylation modification status. We found that in patients with sepsis, multiple metabolic pathways were significantly up-regulated, whereas infiltration by immune cells was significantly decreased. We identified nine potential biomarkers associated with sepsis progression. In addition, we found evidence that the expression of MAPK14, FGR, and CD247 is modified by methylation, which may facilitate the diagnosis of sepsis and septic shock.

Based on analysis of variance, WGCNA and differential expression analysis, we identified genes that may be associated with sepsis and septic shock. Enrichment analysis revealed that the genes involved in inflammatory and metabolic responses were up-regulated in sepsis, whereas the genes corresponding to immune responses were significantly suppressed. This is similar to previous findings.³⁰ Increased galactosylation of the crystallizable immunoglobulin G fragment has been linked to increased complement-dependent cytotoxicity, which in turn may contribute to pediatric sepsis.³¹ In patients with sepsis, increased levels of both C3a and C5a in serum have been reported.³² Monocytes from sepsis patients have significantly higher basal glycolysis than monocytes from healthy controls.³³ Aerobic glycolysis has been found to contribute to sepsis development by regulating inflammasome activation in macrophages.³⁴ Both bacterial and fungal sepsis have been observed to cause a shift in cellular metabolism towards glycolysis, which is associated with an impaired ability of leukocytes to produce pro-inflammatory cytokines upon a second stimulus.³⁵ Moreover, lymphocyte loss, equally affecting B and T cells, was demonstrated in septic patients.^36,37 Accumulating evidence supports that immunosuppression is one of the major contributors to sepsis-related mortality, and that T cell exhaustion is one of the most severe responses.^38,39

Consistent with our KEGG enrichment results, HIF-1α contributes to the pathogenic role of macrophages in sepsis.^40,41 HIF-1α promotes the expression of multiple gene products, including enzymes that promote glycolysis.⁴² Some studies consider HIF as a potential biomarker of sepsis, but this remains controversial.⁴³ Tumor necrosis factor α (TNF-α) contributes to sepsis immunosuppression by increasing apoptosis.⁴⁴ TNF-α levels in plasma increase progressively as sepsis progresses, with the highest values found in patients with septic shock. This finding supports a prognostic biomarker role for this cytokine.⁴⁵ Th17 cytokine production is reduced in sepsis, which may adversely affect long-term mortality.⁴⁶ Treatment with IL-7 increases the responsiveness of Th17 cells and reduces the mortality of secondary fungal infections, making IL-17 a potential therapeutic agent.⁴⁷ A ratio of Th2 to Th1 cells are clearly up-regulated in severe sepsis patients, and their continuous dynamic elevation is associated with acquired infection and 28-day mortality in the intensive care unit.⁴⁸

Patients with sepsis have impaired innate and adaptive antibacterial immunity, which renders them unable to control primary and secondary infections.⁴⁹ Our analysis of immune cell infiltration in sepsis patients showed that the levels of infiltration of some innate immune cells, such as macrophages and mast cells, were significantly higher. The host defends against pathogens by mounting a systemic inflammatory response in innate immune cells, such as macrophages, which produce pro-inflammatory cytokines and chemokines and which initiate coagulation cascade within minutes.⁵⁰ In our clustering of immune cell types based on differential infiltration, we identified four clusters. Each is characterized by unique immune responses that may influence the disease state of the patient.⁵¹

On the other hand, we found 9 key genes that may influence the development of sepsis. Among them, MAPK14, FGR, and CD247 expression may be modified by methylation. MAPK14 has already been identified as a key gene and potential therapeutic target in neonatal sepsis.^52,53 Deficiency of MAPK14 in macrophages protects mice from lipopolysaccharide (LPS)-induced sepsis.⁵⁴ FGR proto-oncogene (FGR) is a Src family kinase that is expressed in innate immune cells, including macrophages and granulocytes, and is considered a major signaling molecule downstream of many immune cell receptors.⁵⁵ FGR ^−/− mice have leukocyte migration disorder, and FGR is involved in the release of proinflammatory mediators.⁵⁶ Studies have shown that CD247 is a potential new biomarker for sepsis and may be beneficial for the prognosis of patients with sepsis.^57,58 CD247 is involved in the immune response to sepsis caused by Staphylococcus aureus infection.⁵⁹

For the other key genes, RHOG has already been linked to sepsis through its participation in inflammatory processes.⁶⁰ The amount of LAT on the cell surface may determine the extent of T cell activation.⁶¹ Systemic inflammation in patients with acute sepsis causes persistent T cell dysfunction, which leads to immunosuppression.³⁷ PRKACB is differentially expressed in neonatal sepsis, and it may influence disease progression via the MAPK signaling pathway.^62,63 ITK has been found to regulate thermal homeostasis during sepsis via its effects on mast cells.⁶⁴ Previous work has found that IL2RB is down-regulated in sepsis, which is consistent with our study,⁶⁵ and that IL2RB is negatively correlates with organ failure and mortality in sepsis.⁶⁶ Our results reveal UBE2Q2 as a potential key gene for sepsis diagnosis, and we are unaware of previous studies showing an association between UBE2Q2 and sepsis.

Our study also has certain limitations. First, the data for this analysis were from public databases and involved small samples, which may lead to biased interpretation of the results. Second, molecular experiments for validation are missing, and follow-up with larger clinical samples is needed to validate our main results. Third, further study is required to confirm whether methylation modifications found in key genes are associated with the progression of sepsis.

Conclusion

The present study demonstrated heightened metabolic response and suppressed immune response in patients with sepsis and septic shock. The nine key genes that we identified may be useful for diagnosing sepsis and monitoring its development, while their specific functions and mechanisms in the disease need further study. Indeed, further experimental analysis is needed to validate our major results.