Bioinformatics Analysis for Multiple Gene Expression Profiles in Sepsis.

BACKGROUND This work aimed to screen key biomarkers related to sepsis progression by bioinformatics analyses. MATERIAL AND METHODS The microarray datasets of blood and neutrophils from patients with sepsis or septic shock were downloaded from Gene Expression Omnibus database. Then, differentially expressed genes (DEGs) from 4 groups (sepsis versus normal blood samples; septic shock versus normal blood samples; sepsis neutrophils versus normal controls and septic shock neutrophils versus controls) were respectively identified followed by functional analyses. Subsequently, protein-protein network was constructed, and key functional sub-modules were extracted. Finally, receiver operating characteristic analysis was conducted to evaluate diagnostic values of key genes. RESULTS There were 2082 DEGs between blood samples of sepsis patients and controls, 2079 DEGs between blood samples of septic shock patients and healthy individuals, 6590 DEGs between neutrophils from sepsis and controls, and 1056 DEGs between neutrophils from septic shock patients and normal controls. Functional analysis showed that numerous DEGs were significantly enriched in ribosome-related pathway, cell cycle, and neutrophil activation involved in immune response. In addition, TRIM25 and MYC acted as hub genes in protein-protein interaction (PPI) analyses of DEGs from microarray datasets of blood samples. Moreover, MYC (AUC=0.912) and TRIM25 (AUC=0.843) had great diagnostic values for discriminating septic shock blood samples and normal controls. RNF4 was a hub gene from PPI analyses based on datasets from neutrophils and RNF4 (AUC=0.909) was capable of distinguishing neutrophil samples from septic shock samples and controls. CONCLUSIONS Our findings identified several key genes and pathways related to sepsis development.

Background

Sepsis represents a serious public health concern and septic shock is a subset of sepsis [1]. More specifically, sepsis is a life-threatening organ dysfunction triggered by a dysregulated host response to infection while septic shock is a deadly health condition caused by uncontrolled sepsis and accompanied by the cellular or metabolic abnormalities [2]. It has been estimated that over 30 million individuals suffer from sepsis worldwide and there has been an increasing incidence and mortality rate during the past decade [3]. Although sepsis is managed with adequate treatment, the overall clinical outcomes have been unsatisfactory. Moreover, existing research suggests that early diagnosis has been considerably difficulty due to several comorbidities and the lack of an effective prediction technique [4]. Therefore, there is a pressing need to identify new biomarkers related to sepsis for the early diagnosis, monitoring, and therapeutic interventions of sepsis.

Analyzing the characterization of patients based on the molecular basis is a powerful approach that can be used for screening promising targets related to disease diagnosis and prognosis. For example, Jekarl et al. [5] evaluated the diagnostic and prognostic values of several cytokines, chemokines, and growth factors. Consequently, their study revealed that hepatic growth factor represented the best advantage for sepsis recognition and epidermal growth factor exhibited a favorable prognosis [5]. Spoto et al. found that the combination of procalcitonin and midregional pro adrenomedullin greatly improved the diagnosis of sepsis [6]. Notably, the bioinformatics analysis based on profiling gene expression data has also been successfully used for identifying the significant biomarkers for various diseases including sepsis [7-9]. Qi et al. screened several differentially expressed genes (DEGs) between sepsis survivors and non-survivors by whole genome profiling and found upregulated type I interferon signaling pathway was closely correlated with the death of sepsis patients [10]. Recently, Qin et al. [11] performed an RNA-sequencing analysis of 3 sepsis patients and 3 healthy volunteers. They noted that several key mRNAs and pathways, such as the T cell receptor signaling pathway and pathways in cancer, might play essential roles in sepsis development [11]. Although existing evidence has identified multiple important gene signatures related to sepsis, the underlying pathogenesis of sepsis has not been fully understood.

In this study, we conducted a comprehensive bioinformatics analysis based on a larger sample size to screen for promising gene markers for sepsis. The microarray datasets for sepsis were retrieved and downloaded from Gene Expression Omnibus (GEO) database. Consequently, we obtained 2 types of gene expression datasets which were respectively generated from blood samples and neutrophil samples. Notably, neutrophils released to blood are recruited to the inflammatory sites to act as significant players for the host defense against microbial invaders [12]. Moreover, studies show that neutrophils exert important roles in infection control in sepsis process and their migratory activity are impaired during sepsis, leading to dysregulated immune responses [13]. Besides, neutrophils also can release many signaling cascades and inflammatory mediators, which will amplify inflammatory responses and eventually cause multiple system organ failure [14]. Sepsis generally tends to develop into septic shock accompanied by vascular dysfunction due to lack of efficient early control for sepsis. Thus, targeting the underlying gene markers in neutrophils will also facilitate identifying key targets for early diagnosis and intervention of sepsis. Herein, for this study, we first identified the DEGs between sepsis or septic shock blood and normal controls, respectively. Then, we respectively conducted DEGs analyses between neutrophil from patients with sepsis or septic shock and normal controls. Furthermore, functional analyses were performed, and the relationships of proteins were also investigated. Finally, the diagnostic values of key genes were evaluated by the receiver operating characteristic (ROC) analyses. This study will identify several biomarkers associated with sepsis in blood and gene markers in neutrophils during sepsis, which will offer deeper insights into the molecular mechanisms of sepsis.

Material and Methods

Data collection

To identify several potential gene signatures associated with sepsis, we performed a systematic bioinformatics analysis based on the gene expression datasets for sepsis, which were retrieved and downloaded from the GEO database () using the searching terms of “Sepsis” AND “Homo sapiens” (porgn) AND “gse” (Filter) [15]. Notably, the included datasets in this study needed to meet the following criteria: 1) the datasets should be the gene expression profiles; and 2) the gene expression data were generated from sepsis cases and normal controls.

Screening DEGs

Firstly, raw gene expression data was pre-processed to eliminate the heterogeneity from different platforms, which primarily included data normalization and log2 transformation. Then, the differential expression analysis between sepsis and normal control groups was carried out. In brief, the P values and false discovery rate (FDR) were respectively calculated using the R metaMA package which has been frequently used for the integrated analysis based on different microarray datasets [16]. Herein, those genes with FDR <0.01 were considered as the DEGs. Furthermore, the bidirectional hierarchical clustering analyses of DEGs extracted were also conducted with the R pheatmap package ().

Functional enrichment analyses

To further recognize the underlying biological functions of identified DEGs in cells, we performed the Gene Ontology (GO) functional analysis to annotate all DEGs according to 3 categories: molecular function (MF), cellular component (CC) and biological process (BP). Meanwhile, the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was also carried out to elucidate DEG pathways. Herein, we used the R clusterProfiler package to undertake GO and KEGG analyses [17]. The parameters selected were as follows: OrgDb=org.Hs.eg.db, pAdjustMethod=“Benjamini & Hochberg”, PvalueCutoff=0.05, qvalueCutoff=0.2.

Protein–protein interaction (PPI) network analysis

Generally, proteins function as team players in a complicated dynamic network. Therefore, identifying and analyzing protein–protein interaction (PPI) pairs promote recognizing key interacting protein partners and better understanding how proteins are organized into specific functional units. Herein, to extract key proteins and explore their potential regulatory roles in the molecular mechanisms of sepsis, we conducted PPI analyses based on the Biological General Repository for Interaction Datasets (BioGRID) database [18]. The Cytoscape software (version 3.0; ) was used to construct the PPI network [19]. The CytoNCA plugin (version 2.1.6; ) was employed to evaluate the topological properties of PPI network nodes according to the parameter of without weight (unweighted) [20]. Accordingly, the scores of degree centrality (DC), betweenness centrality (BC), and closeness centrality (CC) were computed and ranked. Moreover, those nodes with high degree were defined as hub proteins which were functional points and played critical roles in the network structure. Additionally, the MCODE plugin (version 1.5.1; ) of Cytoscape software was used to perform the PPI sub-module analysis using the default parameters of Degree Cutoff=2, Node Score Cutoff=0.2, K-core=2 and Max.Depth=100. Furthermore, the clustering analysis was performed to obtain significantly functional modules and those sub-modules with PPI score >4 were retained. Finally, the KEGG enrichment analysis of DEGs in functional modules was carried out.

ROC analysis

To evaluate whether identified DEGs had important diagnostic values for sepsis or septic shock, the ROC analyses were carried out using the R “pROC” package () [21]. The area under the ROC curve (AUC) of each gene was calculated. The diagnostic accuracy of key biomarkers for sepsis was evaluated with AUC value. Herein, when AUC value was greater than 0.8, the gene could differentiate patients with sepsis or septic shock from healthy individuals.

Results

Acquisition of microarray datasets

We searched the GEO database to obtain the gene expression microarray data about sepsis and found that these profiles were generated by blood samples or neutrophil samples from patients with sepsis or septic shock. In this study, we first identified the potential biomarkers related to sepsis or septic shock by analyzing the gene expression profiles from blood samples. Seven datasets (GSE69528, GSE46955, GSE54514, GSE32707, GSE28750, GSE13015, and GSE9960) from blood samples of 341 sepsis patients and 145 normal controls were considered as group 1 and used to identify underlying gene markers involved with sepsis (Table 1). Three datasets obtained from blood samples of 204 patients with septic shock and 89 healthy individuals, including GSE33118, GSE95233, and GSE57065, were regarded as the group 2 to extract key biomarkers correlated with septic shock (Table 1). Notably, overwhelming evidence has demonstrated that neutrophils play pivotal roles in the pathogenesis of sepsis. Therefore, we next focused on targeting neutrophils to screen the potential gene signatures associated with sepsis or septic shock, which will facilitate developing effective medical treatment strategies for conquering sepsis or controlling its development. Accordingly, 3 datasets (GSE49755, GSE49756, and GSE49757) were used as group 3, involving neutrophil samples from 88 sepsis and 48 healthy controls (Table 1). Meanwhile, GSE123729 and GSE64457 as group 4 were generated based on neutrophil samples from 30 septic shock patients and 19 normal controls (Table 1).

Table 1

Gene expression datasets for sepsis and septic shock.

GEO accession	Type	Case	Control	Platform	Author
GSE69528	Blood samples (sepsis vs. controls)	83	24	GPL10558 IlluminaHumanHT-12 V4.0 expression beadchip	Scott Presnell
GSE46955	Blood samples (sepsis vs. controls)	8	6	GPL6104 IlluminahumanRef-8 v2.0 expression beadchip	Michael Poidinger
GSE54514	Blood samples (sepsis vs. controls)	127	36	GPL6947 IlluminaHumanHT-12 V3.0 expression beadchip	Grant Parnell
GSE32707	Blood samples (sepsis vs. controls)	30	34	GPL10558 IlluminaHumanHT-12 V4.0 expression beadchip	Judie Ann Howrylak
GSE28750	Blood samples (sepsis vs. controls)	10	20	GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array	Gareth Price
GSE13015	Blood samples (sepsis vs. controls)	29	5	GPL6947 IlluminaHumanHT-12 V3.0 expression beadchip	Damien Chaussabel
GSE9960	Blood samples (sepsis vs. controls)	54	16	GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array	Benjamin Man Piu Tang
GSE33118	Blood samples (septic shock vs. controls)	20	42	GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array	Wolfgang Raffelsberger
GSE95233	Blood samples (septic shock vs. controls)	102	22	GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array	Julien Textoris
GSE57065	Blood samples (septic shock vs. controls)	82	25	GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array	Marie-Angélique Cazalis
GSE49755	Neutrophil samples (sepsis vs. controls)	29	17	GPL10558 IlluminaHumanHT-12 V4.0 expression beadchip	Damien Chaussabel
GSE49756	Neutrophil samples (sepsis v controls)	24	12	GPL10558 IlluminaHumanHT-12 V4.0 expression beadchip	Damien Chaussabel
GSE49757	Neutrophil samples (sepsis vs. controls)	35	19	GPL10558 IlluminaHumanHT-12 V4.0 expression beadchip	Damien Chaussabel
GSE123729	Neutrophil samples (septic shock vs. controls)	15	11	GPL21970 [HuGene-2_0-st] Affymetrix Human Gene 2.0 ST Array	Carsten Sticht
GSE64457	Neutrophil samples (septic shock vs. controls)	15	8	GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array	Julien Textoris

GEO – Gene Expression Omnibus.

Identification of DEGs in sepsis and septic shock

To eliminate the heterogeneity from different sequencing platforms, the R metaMA package was used to identify the DEGs in different groups using the FDR <0.01 as the screening criterion. Here, we extracted 2082 DEGs from group 1, including 928 upregulated genes and 1104 downregulated genes. There were 2079 DEGs (865 upregulated genes and 1214 downregulated genes) from group 2. For group 3, 6590 DEGs were obtained, which contained 3070 upregulated genes and 3520 downregulated genes. We also identified 1056 DEGs (478 upregulated genes and 578 downregulated genes) from group 4. Furthermore, the bidirectional clustering analysis of the DEGs in different groups suggested that they remarkably separated patients with sepsis or septic shock from healthy controls, which further implied that these DEGs may be closely associated with sepsis (Figure 1).

Figure 1

The bidirectional clustering analysis of differentially expressed genes (DEGs). (A) The heatmap of DEGs from blood samples of sepsis and healthy controls. (B) The heatmap of DEGs from blood samples of septic shock and healthy controls. (C) The heatmap of DEGs from neutrophil samples of sepsis and normal controls. (D) The heatmap of DEGs from neutrophil samples of septic shock and normal controls. The “label” shows the type of microarray dataset. Each row and column respectively represent the gene and sample. The blue-green color represents the control samples and the pink color represents the sepsis or septic shock samples.

We further performed the functional analyses for DEGs in different groups to explore their underlying biological roles. The DEGs in group 1 were enriched in 266 GO-BP terms and several significant terms were exhibited in Figure 2A such as ribosome biogenesis, ribonucleoprotein complex biogenesis, and ncRNA processing. Meanwhile, there were 97 GO-CC terms; specific granule, secretory granule lumen, and tertiary granule were top 3 remarkable GO-CC terms (Figure 2A). The GO-MF analysis also revealed that these DEGs in group 1 were dramatically enriched in 10 key GO-MF terms such as rRNA binding, catalytic activity, acting on RNA and translation factor activity, RNA binding (Figure 2A). In addition, they also were responsible for 15 KEGG pathways, which mainly included ribosome, ribosome biogenesis in eukaryotes and viral carcinogenesis (Figure 3A).

Figure 2

The GO functional annotation of DEGs. (A) The top 10 GO terms enriched by DEGs from blood samples of sepsis and healthy controls in 3 GO categories (GO-BP; GO-CC and GO-MF). (B) The top 10 GO terms enriched by DEGs from blood samples of sepsis and healthy controls in 3 GO categories (GO-BP; GO-CC and GO-MF). (C) The top 10 GO terms enriched by DEGs from neutrophil samples of sepsis and normal controls in 3 GO categories (GO-BP; GO-CC and GO-MF). (D) The top 10 GO terms enriched by DEGs from neutrophil samples of septic shock and normal controls in 3 GO categories. GO – Gene Ontology; DEGs – differentially expressed genes; MF – molecular function; CC – cellular component; BP – biological process.

Figure 3

The KEGG pathway enrichment analysis of DEGs. (A) The top 15 KEGG pathways of DEGs from blood samples of sepsis patients and healthy controls. (B) The top 15 KEGG pathways of DEGs from blood samples of sepsis patients and healthy controls. (C) The top 15 KEGG pathways of DEGs from neutrophil samples of sepsis and normal controls. (D) The significant enriched KEGG pathway of DEGs from neutrophil samples of septic shock and normal controls. The color close to red shows a higher significance. KEGG – Kyoto Encyclopedia of Genes and Genomes; DEGs – differentially expressed genes.

We found that the DEGs in group 2 were associated with 585 GO-BP terms. Moreover, multiple GO-BP terms associated with neutrophil activity were markedly enriched by numerous DEGs, including neutrophil activation, neutrophil activation involved in immune response, neutrophil degranulation, and neutrophil mediated immunity (Figure 2B). Additionally, these DEGs played essential roles in 101 GO-CC terms containing 3 most enriched terms: specific granule, secretory granule lumen, and cytoplasmic vesicle lumen; Figure 2B). We also observed that these DEGs were involved in 22 GO-MF terms such as catalytic activity acting on RNA, catalytic activity acting on a tRNA and ligase activity (Figure 2B). The KEGG analysis suggested that DEGs in group 2 were also significantly enriched in 15 KEGG pathways, such as Th17 cell differentiation, Th1 and Th2 cell differentiation, and primary immunodeficiency (Figure 3B).

For DEGs extracted from group 3, the GO-BP analysis showed that they were involved in 354 terms and several neutrophil activity-related terms such as neutrophil activation, neutrophil activation involved in immune response, neutrophil degranulation and neutrophil mediated immunity were most remarkably enriched (Figure 2C). These DEGs were also strongly correlated with 83 GO-CC terms, primarily including secretory granule membrane, specific granule, and specific granule membrane (Figure 2C). Meanwhile, GO-MF analysis of DEGs in group 3 indicated that they exerted crucial roles in 15 protein activity and binding-related terms such as protein serine/threonine kinase activity, 14-3-3 protein binding, phosphatase binding, and Rho GTPase binding (Figure 2C). Besides, these DEGs were enriched in 19 KEGG pathways, mainly including insulin signaling pathway, Fc gamma R-mediated phagocytosis, and platelet activation (Figure 3C).

For those DEGs in group 4, we noted that they were closely related to 242 GO-BP terms and significantly enriched terms contained multiple neutrophil activity-associated terms such as neutrophil activation, neutrophil degranulation, neutrophil mediated immunity and neutrophil activation involved in immune response (Figure 2D). Our functional annotation analyses also demonstrated that these DEGs were enriched in 63 GO-CC terms and 11 GO-MF terms. As displayed in Figure 2D, tertiary granule, ficolin-1-rich granule, and nuclear speck were 3 remarkable terms while protein binding-related terms such as histone binding, ubiquitin-like protein ligase binding, and ubiquitin protein ligase binding were dramatically enriched by DEGs from group 4 (Figure 2D). Moreover, the KEGG analysis revealed that there was only one KEGG pathway (autophagy-animal pathway; Figure 3D).

PPI analyses

PPI network analysis of DEGs in group 1 revealed that there were 1600 nodes and 8911 interaction pairs (Supplementary Figure 1). In this study, the top 5 genes according to the DC value were regarded as the hub genes. Herein, the hub genes in the PPI network of DEGs from group 1 contained upregulated tripartite motif containing 25 (TRIM25) and amyloid beta precursor pro (APP), and downregulated ring finger protein 4 (RFN4), MYC proto-oncogene (MYC), and exportin 1 (XPO1) as showed in Supplementary Table 1. Then, the functional sub-modules were further extracted. Our results suggested that there were 4 significant sub-modules (Figure 4A–4D). Specifically, there were 34 genes (30 upregulated genes and 4 downregulated genes) in sub-module 1, primarily containing several ribosomal protein (RP) genes such as RPL3, RPL8, and RPL12. We also found that TRIM25 and MYC were key genes in sub-module 1. The functional enrichment analysis showed that several RP genes, such as RPL3, RPL8, RPL12, and RPS18, were significantly involved in ribosome pathway while MYC participated in the cancer-related pathways such as thyroid cancer pathway (Supplementary Table 2). For sub-module 2, there were 15 downregulated genes including eukaryotic translation initiation factor 3 (EIF3) complex gene (EIF3A-B, EIF3D, EIF3F-I, EIF3K-M) and eukaryotic translation initiation factor 4B (EIF4B), RNF4, CD81 molecule (CD81), poly(A) binding protein cytoplasmic 1 (PABPC1), and HECT and RLD domain containing E3 ubiquitin protein ligase 2 (HERC2). Our KEGG analysis revealed that numerous EIF genes (EIF3A-B, EIF3D, EIF3F-I, and EIF4B) primarily participated in RNA transport pathway (Supplementary Table 2). Moreover, PABPC1 was strongly associated with RNA pathway including RNA transport, RNA degradation, and mRNA surveillance pathway (Supplementary Table 2). There were 29 genes (8 upregulated genes and 21 downregulated genes) in sub-module 3. The KEGG enrichment analysis showed that several genes were predominately enriched in RNA degradation (exosome component/EXOSC 1, namely EXOSC1, and EXOSC4-10), spliceosome (eukaryotic translation initiation factor 4A3/EIF4A3, pre-mRNA processing factor 19/PRPF19, splicing factor 3A subunit 1/SF3A1, SF3A2, and SF3A3), and ribosome biogenesis in eukaryotes (uridine triphosphate 14 small subunit processome component/UTP14A and GTP binding protein 4/GTPBP4; Supplementary Table 2). Besides, we found 13 genes (3 upregulated genes and 10 downregulated genes) in PPI sub-module 4 and top 2 KEGG pathways were respectively spliceosome and RNA degradation pathways (Supplementary Table 2). Furthermore, the LSm family of RNA-binding proteins LSM4 and LSM8 played pivotal roles in these 2 pathways (Supplementary Table 2).

Figure 4

The functional sub-network analysis of PPI network based on the blood and neutrophil samples from sepsis patients. (A–D) The 4 sub-modules from PPI network of DEGs from blood samples of sepsis and healthy controls. (E, F) The 3 sub-modules from PPI network of DEGs from neutrophil samples of sepsis and healthy controls. The red color nodes represent upregulated genes and the blue nodes show the downregulated genes. PPI – protein–protein interaction; DEGs – 0differentially expressed genes.

Our results indicated 5145 gene nodes and 46251 protein pairs in PPI network of DEGs in group 2 (Supplementary Figure 2). Downregulated KIAA1429/VIRMA (vir like M6A methyltransferase associated), upregulated TRIM25, upregulated RecQ-like helicase 4 (RECQL4), downregulated MYC, and downregulated ELAV-like RNA binding protein 1 (ELAVL1) were considered as the key hub genes (Supplementary Table 1). Furthermore, there were 8 PPI sub-modules (Figure 5A–5H). The sub-module 1 contained 34 genes (3 upregulated genes and 31 downregulated genes). More specifically, numerous RP genes such as RPL3, RPL8, and RPS3 were clustered in sub-module 1. However, no significant KEGG pathway was found (Supplementary Table 2). The sub-module 2 consisted of 28 upregulated genes and 26 downregulated genes. These genes were mainly associated with cell cycle (tumor protein P53/TP53 and cell division cycle 25B/CDC25B), MAPK signaling pathway (TP53, CDC25B, MYC associated factor X/MAX and arrestin beta 2/ARRB2), and microRNAs (miRNAs) in cancer pathways (TP53, CDC25B, histone deacetylase 4/HDAC4, and heterogeneous nuclear ribonucleoprotein K/HNRNPK; Supplementary Table 2). There were 3 upregulated genes and 23 downregulated genes in PPI sub-module 3. Moreover, the KEGG enrichment analysis suggested that spliceosome, oocyte meiosis, and cell cycle were significantly enriched pathways. Seven genes (heterogeneous nuclear ribonucleoprotein A3/HNRNPA3, U2 small nuclear RNA auxiliary factor 1/U2AF1, SF3A1, pre-mRNA processing factor 19/PRPF19, small nuclear ribonucleoprotein polypeptide A/SNRPA1, DEAH-box helicase 15/DHX15, and DEAD-box helicase 5/DDX5) were involved in the spliceosome pathway. Cell division cycle 25C (CDC25C) and tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein theta (YWHAQ) were enriched in oocyte meiosis and cell cycle pathways (Supplementary Table 2). Four upregulated genes and 7 downregulated genes were found in sub-module 4. Three pathways (spliceosome, ribosome, and PI3K-Akt signaling pathway) were markedly enriched. Two genes poly (RC) binding protein 1/PCBP1 and heterogeneous nuclear ribonucleoprotein/HNRNP were implicated with spliceosome pathway. RPS27L and RPS28 had key roles in ribosome pathway while integrin subunit alpha 4/ITGA4 and breast cancer 1/BRCA1 were dramatically enriched in PI3K-Akt signaling pathway (Supplementary Table 2). We found that 65 genes (18 upregulated genes and 47 downregulated genes) in sub-module 5 and several genes were significantly enriched in cell cycle (BUB3 mitotic checkpoint protein/BUB3), RNA degradation (EXOSC4-8 and DIS3 like exosome 3’-5’ exoribonuclease/DIS3L), and spliceosome pathway (SF3A2, serine and arginine rich splicing factor 7/SRSF7, PRPF3, small nuclear ribonucleoprotein U5 subunit 40/SNRNP40, peptidylprolyl isomerase E/PPIE and small nuclear ribonucleoprotein polypeptide E/SNRPE; Supplementary Table 2). The sub-module 6 included 13 upregulated genes and 22 downregulated genes were in sub-module 6. The KEGG pathway analysis showed that these genes were mainly responsible for cell cycle, endocrine resistance, and PI3K-Akt signaling pathway. Notably, cyclin dependent kinase 4 (CDK4) was closely associated with these 3 pathways (Supplementary Table 2). For sub-module 7, there were 99 genes, including 19 upregulated genes and 80 downregulated genes. Functional analysis revealed that they predominately played important roles in DNA replication, RNA transport and proteasome. EIF1B, EIF3A, EIF3C, EIF3D, and EIF3J were enriched in RNA transport. Meanwhile proteasome 26S subunit, ATPase 4 (PSMC4), proteasome subunit beta 5 (PSMB5), proteasome activator subunit 1 (PSME1), and proteasome 26S subunit, non-ATPase 1 (PSMD1) were involved in proteasome (Supplementary Table 2). There were 41 upregulated genes and 73 downregulated genes in sub-module 8 and they were mainly enriched metabolism-related pathways such as purine and pyrimidine metabolism. RNA polymerase I and III subunit D (POLR1D), RNA polymerase III subunit C (POLR3C) and POLR3D participated in these metabolic pathways.

Figure 5

The functional sub-network analysis of PPI network based on the blood and neutrophil samples from septic shock patients. (A–H) The 8 sub-modules from PPI network of DEGs from blood samples of septic shock and healthy controls. (I) The one sub-module from PPI network of DEGs from neutrophil samples of septic shock and normal controls. The red color nodes represent upregulated genes and the blue nodes show the downregulated genes. PPI – protein–protein interaction; DEGs – differentially expressed genes.

For the PPI analysis of DEGs in group 3, we found 1375 nodes and 6198 interactions in PPI network (Supplementary Figure 3). Moreover, downregulated COP9 signalosome subunit 5 (COPS5), downregulated methylphosphate capping enzyme (MEPCE), downregulated histone deacetylase 1 (HDAC1), upregulated TRIM25, and downregulated RNF4 acted as the critical hub genes (Supplementary Table 1). Additionally, 3 sub-modules were identified from PPI network (Figure 4E–4G). As indicated in Supplementary Table 2, the sub-module 1 contained 7 upregulated genes and 4 downregulated genes. Many genes were markedly enriched in transcriptional mis-regulation in cancer (B-cell lymphoma 6/BCL6 and zinc finger and BTB domain containing 16/ZBTB16), nicotinate and nicotinamide metabolism (nicotinamide adenine dinucleotide kinase/NADK), acute myeloid leukemia (ZBTB16), Fc epsilon RI signaling pathway (growth factor receptor bound protein 2-associated binding protein 2 (GAB2) and chronic myeloid leukemia (GAB2). There were 11 genes in sub-module 2. Moreover, they were primarily related to endocytosis, protein processing in endoplasmic reticulum pathway, endocrine and other factor-regulated calcium reabsorption, legionellosis, and synaptic vesicle cycle. We noted that clathrin light chain B (CLTB) and adaptor related protein complex 2 subunit Mu 1 (AP2M1) were strongly related to 3 pathways (endocytosis, endocrine and other factor-regulated calcium reabsorption and synaptic vesicle cycle (Supplementary Table 2). Finally, 3 upregulated genes and 3 downregulated genes were clustered in sub-module 3. In addition, the hippo signaling pathway and multiple cancer-associated pathways were remarkably enriched. Ras association domain family member 1 (RASSF1) played critical roles in these significant pathways.

The PPI analysis of DEGs in group 4 revealed that there were 791 nodes and 2548 protein pairs (Supplementary Figure 4). Upregulated ubiquity in C (UBC), downregulated ring finger protein 2 (RNF2), downregulated RNF4, downregulated SUZ12 poly comb repressive complex 2 subunit (SUZ12), and downregulated CYLD lysine 63 deubiquitinase (CYLD) were regarded as hub genes (Supplementary Table 1). Moreover, only one sub-module was extracted, which included 5 downregulated genes (mediator complex subunit 28/MED28, MED29, MED30, MED31, and cyclin dependent kinase 19/CDK19; Figure 5I). Similarly, one significant KEGG pathway (thyroid hormone signaling pathway) was obtained (Supplementary Table 2).

The ROC curve analysis was carried out to assess the diagnostic values of several key DEGs. We found that TRIM25 and MYC was respectively upregulated and downregulated. Furthermore, they also served as the overlapping hub genes in PPI networks of DEGs from group 1 and 2. In addition, RNF4 was downregulated in 4 groups and also an overlapping hub gene in PPI networks of DEGs from group 3 and 4. Herein, we assessed the diagnostic implication of these 3 genes in sepsis. The results indicated that MYC (AUC=0.912) could dramatically differentiate septic shock blood samples and normal control samples. However, this gene exhibited poor identification capability for sepsis blood samples and corresponding controls (AUC=0.712), which suggested that they had important diagnostic values for the detection of septic shock (Figure 6A, 6B). Similarly, our finding revealed that TRIM25 (AUC=0.843) had significant diagnostic value for discriminating septic shock blood samples and normal controls (Figure 6C). Additionally, RNF4 (AUC=0.909) was capable of distinguishing neutrophil samples from septic shock samples and controls (Figure 6D).

Figure 6

ROC curve of 3 differentially expressed genes. The ROC curves to show the diagnostic value of 3 key genes (MYC, TRIM25, and RNF4) in sepsis or septic shock with sensitivity and 1-specificity. The x-axis shows 1-specificity and y-axis represents sensitivity. (A) ROC curve of MYC in sepsis. (B) ROC curve of MYC in septic shock. (C) ROC curve of TRIM25 in septic shock. (D) ROC curve of RNF4 in septic shock. The gene has a great diagnostic value for septic shock when AUC value of this gene >0.8. ROC – receiver operating characteristic; AUC – area under the ROC curve.

Discussion

In this study, we performed a bioinformatics analysis based on the relevant microarray datasets about sepsis to identify several potential markers associated with sepsis. Our findings suggested that ribosome-related pathway, cell cycle, and neutrophil activation involved in immune response might be critical for the sepsis progression. Additionally, TRIM25, RNF4, and MYC were key hub genes in PPI analysis. Moreover, MYC and TRIM25 could significantly discriminate septic shock blood samples and normal control samples. RNF4 had great diagnostic value for differentiating neutrophil samples from septic shock samples and controls.

The functional analysis of DEGs in group 1 revealed that ribosomal biological function played essential roles in the development of sepsis. Furthermore, the KEGG analysis of DEGs in PPI sub-network also showed that ribosome-related pathway was one of significantly enriched pathways. Ribosomes are remarkably crucial for catalyzing protein synthesis in the course of all life forms. Functional ribosomes are ribonucleoprotein complexes assembled by rRNA and proteins. Therefore, the ribosomal dysfunction can lead to various human diseases [22,23]. An early research indicated that sepsis dramatically decreased the formation of 40S ribosome initiation complex in skeletal muscle, which would be not conducive to protein synthesis [24]. Orellana et al. pointed out that sepsis could retard the muscle protein synthesis (MPS) in neonatal pigs via diverse ribosomal molecular mechanisms. For example, the MPS reduction triggered by lipopolysaccharide was closely correlated with the decreased ribosomal efficiency during sepsis [25]. Additionally, our results showed that the cell cycle regulation also might be involved in the progression of sepsis. Similarly, Real et al. examined the expression of RNA transcripts (mRNAs and miRNAs) in exosomes of patients undergoing septic shock and found that exosomes conveyed miRNAs and mRNAs associated with cell cycle regulation [26]. A later study suggested that ghrelin could increase the expression of cell cycle positive regulators and decrease the expression of a cell cycle negative regulator in mice with sepsis, thereby inducing the CD4 T cell proliferation [27]. Our functional analysis also suggested that neutrophil activation pathway was possibly implicated with the sepsis development. Overwhelming evidence has implied that neutrophil activation was significant for sepsis detection [28]. Törnblom et al. argued that increased plasma activin A denoted the activation of neutrophil and neutrophil was accumulated into impaired kidneys, which further indicating neutrophil activation had pivotal pathophysiological roles for septic patients with acute kidney injury [29]. Taken together, further focusing on the ribosome-related pathway, cell cycle and neutrophil activation pathway would contribute to developing promising strategy for conquering sepsis.

It is reported that MYC as a proto-oncogene exerted several roles in cellular processes, such as cell cycle progression, cell proliferation and apoptosis. Liu et al. previously performed a bioinformatics analysis using the GSE26440 dataset, which was obtained from children suffering from septic shock and healthy controls. As a result, they found that MYC was upregulated and NFκB was probably implicated with the progression of septic shock via elevating the expression levels of MYC and other key genes, which was supported by a later study that found that MYC was upregulated and was possibly involved in the molecular mechanism of sepsis [30,31]. However, our finding indicated that MYC was downregulated and it acted as key hub gene in PPI network analysis. A plausible explanation for expression difference of this gene was the different sample source. Zhang et al. analyzed long non-coding RNA (lncRNA) and mRNA sequencing data from rat models of septic shock-induced myocardial depression and they highlighted that lncRNA rPvt1 knockdown elevated the levels of c-Myc, Bax and caspase-3, promoting cell apoptosis in heart tissue of septic rats, suggesting that MYC may be also correlated with septic shock-induced complication [32]. Interestingly, we also observed that MYC exhibited a good discriminatory power for septic shock blood samples and controls. However, this gene could not distinguish sepsis blood samples and normal controls primarily due to limited number of sample size. Taken together, these evidences suggested that MYC was strongly related to the development of sepsis. Herein, our results also demonstrated that TRIM25 was a key gene for sepsis according to the transcriptome profiling from blood samples of sepsis patients. TRIM25 is a key member of tripartite motifs of E3 ubiquitin ligase enzyme which was implicated with numerous biological processes such as the regulation of innate immune responses against viruses [33]. For example, Liu et al. recently revealed that the activation of TRIM25-mediated retinoic acid-inducible gene I in macrophages increased the antiviral immune processes [34]. Moreover, extensive studies have reported that sepsis acted as a proper balance between competing pro- and anti-inflammatory pathways [35,36]. However, there is no evidence that TRIM25 participate in the pathogenesis of sepsis by regulating immune responses. In addition, our ROC analysis showed that TRIM25 had good diagnostic values for predicting septic shock blood samples from normal controls. Therefore, we inferred that MYC and TRIM25 were essential gene markers for the recognizing sepsis and predicting septic shock.

By microarray data analysis from neutrophils, we found that RNF4 was an overlapping key hub gene in PPI network of DEGs between 2 groups (sepsis versus normal controls and septic shock versus normal controls). The ROC analysis suggested that this gene could distinguish septic shock neutrophil samples from normal controls. RNF4 is a transcriptional cofactor containing small ubiquitin-like modifiers (SUMOs)-interacting motif and a RING finger domain which functions as an E3 ubiquitin ligase to link SUMOylation to ubiquitination [37]. Interestingly, a previous study showed that silencing SUMO-targeted ubiquitin ligase RNF4 could increase the levels of SUMOylated c-Myc [38]. However, no report investigated the potential interactive roles of RNF4-MYC-TRIM25 in the molecular mechanism of sepsis.

There are still limitations in the present work. First, the relevant experimental assays should be conducted to confirm the underlying biological roles of key genes and pathways in sepsis. Second, the bioinformatics analysis based on the regulatory network was needed to screen more potential gene biomarkers involved in sepsis. Third, the corresponding clinical information also needed to be integrated into further analysis.

Conclusions

Three key pathways (ribosome-related pathway, cell cycle, and neutrophil activation involved in immune response) were closely associated with the pathogenesis of sepsis. Moreover, TRIM25, RNF4, and MYC were key biomarkers related to sepsis progression and acted as significant predicators for septic shock. However, further experimental analysis was still required to validate these findings.

PPI network of DEGs between sepsis patients and healthy controls according the gene expression profile from blood samples. The red color nodes represent upregulated genes and the blue nodes show the downregulated genes. PPI – protein–protein interaction; DEGs – differentially expressed genes.

PPI network of DEGs between septic shock patients and healthy controls according the gene expression profile from blood samples. The red color nodes represent upregulated genes and the blue nodes show the downregulated genes. PPI – protein–protein interaction; DEGs – differentially expressed genes.

PPI network of DEGs between sepsis patients and healthy controls according the gene expression profile from neutrophil samples. The red color nodes represent upregulated genes and the blue nodes show the downregulated genes. PPI – protein–protein interaction; DEGs – differentially expressed genes.

PPI network of DEGs between septic shock patients and healthy controls according the gene expression profile from neutrophil samples. The red color nodes represent upregulated genes and the blue nodes show the downregulated genes. PPI – protein–protein interaction; DEGs – differentially expressed genes.

Supplemenatry Table 1

The hub genes in the PPI network.

Group 1			Group 2			Group 3			Group 4
Gene	Degree	Type	Gene	Degree	Type	Gene	Degree	Type	Gene	Degree	Type
TRIM25	329	Up	KIAA1429	1056	Down	TRIM25	313	Up	RNF4	112	Down
MYC	317	Down	TRIM25	796	Up	RNF4	185	Down	UBC	111	Up
APP	269	Up	MYC	713	Down	COPS5	110	Down	RNF2	77	Down
XPO1	206	Down	ELAVL1	697	Down	MEPCE	91	Down	SUZ12	66	Down
RNF4	203	Down	RECQL4	520	Up	HDAC1	89	Down	CYLD	60	Down

Group 1: the PPI network analysis of DEGs between sepsis and normal controls from blood samples; Group 2: the PPI network analysis of DEGs between septic shock and normal controls from blood samples; Group 3: the PPI network analysis of DEGs between sepsis and normal controls from neutrophil samples; Group 4: the PPI network analysis of DEGs between septic shock and normal controls from neutrophil samples. PPI – protein–protein interaction; DEGs – differentially expressed genes; Up represents the upregulation of the gene and Down represents the downregulation of the gene.

Supplementary Table 2

The Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of differentially expressed genes (DEGs) in protein–protein interaction (PPI) sub-modules (top 5).

Names	Terms	Gene count	P-value	Gene symbols
Group1_sub-module 1	Ribosome	28	2.78E-62	PRL19, RPL37A, RPL23A, RPL22, RPL3, RPL32, RPL35A, RPL7, PRS18, RPS19, RPS12, RPS23, RPS3, RPS15A, RPL30, RPL14, RPL38, RPL15, RPL10A, RPL8, RPL12, RPL11, RPS16, RPS4X, RPS5, RPS28, RPS25
	Thyroid cancer	1	0.025	MYC
	Bladder cancer	1	0.035	MYC
	Endometrial cancer	1	0.044	MYC
	Acute myeloid leukemia	1	0.048	MYC
Group1_sub-module 2	RNA transport	9	3.35E-18	PABPC1, EIF3A, EIF3B, EIF3D, EIF3F, EIF3G, EIF3H, EIF3I, EIF4B
	Malaria	1	0.0187	CD81
	B cell receptor signaling pathway	1	0.0276	CD81
	RNA degradation	1	0.0290	PABPC1
	mRNA surveillance pathway	1	0.0345	PABPC1
Group1_sub-module 3	RNA degradation	9	6.46E-18	EXOSC1, EXOSC4, EXOSC5, EXOSC6, EXOSC7, EXOSC8, EXOSC9, EXOSC10, SKIV2L2
	Spliceosome	5	5.37E-08	EIF4A3, PRPF19, SF3A1, SF3A2, SF3A3
	Ribosome biogenesis in eukaryotes	2	0.0020	UTP14A, GTPBP4
	Rap1 signaling pathway	2	0.0105	TIAM1, SIPA1L1
	MicroRNAs in cancer	2	0.0202	CDC25B, HDAC4
Group1_sub-module 4	Spliceosome	3	1.11E-05	LSM4, LSM8, PPIH
	RNA degradation	2	0.0003	LSM4, LSM8
	RNA polymerase	1	0.0107	POLR2D
	Basal transcription factors	1	0.0149	TAF1
	Pyrimidine metabolism	1	0.0341	POLR2D
Group2_sub-module 1	–	–	–
Group2_sub-module 2	Cell cycle	2	0.0128	TP53, CDC25B
	MAPK signaling pathway	4	0.0003	TP53, CDC25B, MAX, ARRB2
	MicroRNAs in cancer	4	0.0007	TP53, CDC25B, HDAC4, HNRNPK
	Viral carcinogenesis	3	0.0029	TP53, HDAC4, HNRNPK
	Chronic myeloid leukemia	2	0.0047	TP53, GAB2
Group2_sub-module 3	Spliceosome	7	3.76E-12	HNRNPA3, U2AF1, SF3A1, PRPF19, SNRPA1, DHX15, DDX5
	Oocyte meiosis	2	0.0030	CDC25C, YWHAQ
	Cell cycle	2	0.0031	CDC25C, YWHAQ
	Pathogenic Escherichia coli infection	1	0.0360	YWHAQ
	Shigellosis	1	0.0423	U2AF1
Group2_sub-module 4	Spliceosome	2	0.0006	PCBP1, HNRNP
	Ribosome	2	0.0007	RPS27L, RPS28
	PI3K-Akt signaling pathway	2	0.0039	ITGA4, BRCA1
	Intestinal immune network for IgA production	1	0.0140	ITGA4
	Legionellosis	1	0.0154	EEF1A1
Group2_sub-module 5	Cell cycle	1	0.1849	BUB3
	RNA degradation	6	5.10E-09	EXOSC6, DIS3L, EXOSC4, EXOSC8, EXOSC5, EXOSC7
	Spliceosome	6	1.18E-07	SRSF7, PRPF3, SNRNP40, PPIE, SNRPE, SF3A2
	Proteasome	4	1.19E-06	PSMD6, PSMD7, PSMB6, PSMD12
	Ribosome	5	3.86E-06	MRPL10, MRPS14, MRPL30, MRPL4, MRPL28
Group2_sub-module 6	Axon guidance	3	0.0005	CFL1, KRAS, EPHA2
	Cell cycle	2	0.0055	CDKN2C, CDK4
	Endocrine resistance	3	9.50E-05	CDKN2C, CDK4, KRAS
	PI3K-Akt signaling pathway	4	0.0002	KRAS, CDK4, FGFR4, EPHA2
	Signaling pathways regulating pluripotency of stem cells	3	0.0003C	SOX2, KRAS, FGFR4
Group2_sub-module 7	DNA replication	3	0.0001	RFC5, RFC1, SSBP1
	RNA transport	6	5.79E-06	EIF1B, EIF3A, EIF3C, EIF3D, EIF3J, SRRM1
	Proteasome	4	6.41E-06	PSMC4, PSMB5, PSME1, PSMD1
	Nucleotide excision repair	4	8.18E-06	RFC5, RFC1, XPC, ERCC6
	Cell cycle	1	0.2675	CCND1
Group2_sub-module 8	Cell cycle	4	0.0005	E2F1, SMC1A, TFDP1, MCM3
	Metabolic pathways	19	3.30E-09	DNMT3A, CYC1, B3GNT2, ALDH18A1, PFAS, SORD, MOGS, SMS, POLR3D, MAN2A1, NAGK, POLE3, POLR1D, LDHAL6B, POLR3C, SDHA, NTPCR, PFKP, DLST
	Purine metabolism	7	1.04E-06	PNPT1, PFAS, POLR3D, NTPCR, POLE3, POLR1D, POLR3C
	Pyrimidine metabolism	5	1.68E-05	POLE3, POLR3C, PNPT1, POLR1D, POLR3D
	RNA polymerase	3	0.0001	POLR3C, POLR1D, POLR3D
Group3_sub-module 1	Transcriptional misregulation in cancer	2	0.0011	BCL6, ZBTB16
	Nicotinate and nicotinamide metabolism	1	0.0085	NADK
	Acute myeloid leukemia	1	0.0159	ZBTB16
	Fc epsilon RI signaling pathway	1	0.0189	GAB2
	Chronic myeloid leukemia	1	0.0203	GAB2
Group3_sub-module 2	Endocytosis	4	6.04E-07	CLTB,ARF1,IGF2R,AP2M1
	Protein processing in endoplasmic reticulum	3	1.21E-05	SEC24C,SEC23B,MAN1A1
	Endocrine and other factor-regulated calcium reabsorption	2	8.12E-05	CLTB,AP2M1
	Legionellosis	2	0.0001	RAB1A,ARF1
	Synaptic vesicle cycle	2	0.0001	CLTB,AP2M1
Group3_sub-module 3	Hippo signaling pathway -multiple species	3	9.45E-09	SAV1, RASSF1,LATS2
	Hippo signaling pathway	3	1.20E-06	SAV1, RASSF1,LATS2
	Ras signaling pathway	2	0.0004	RASSF1,MET
	MicroRNAs in cancer	2	0.0008	RASSF1,MET
	Pathways in cancer	2	0.0015	RASSF1,MET
Group4_sub-module	Thyroid hormone signaling pathway	1	0.0149	MED30

Group 1: the blood samples of sepsis and healthy controls; Group 2: the blood samples of septic shock and healthy controls; Group 3: the neutrophil samples of sepsis and normal controls; Group 4: the neutrophil samples of septic and normal controls.