Cases in the cohort

This was a retrospective study, and our protocol was approved by the Institutional Review Board (IRB) of Incheon St. Mary’s Hospital and Seoul St. Mary’s Hospital in accordance with the Declaration of Helsinki. As public datasets were used, informed consent was waived by the IRB. The public datasets were downloaded from the Gene Express Omnibus (GEO) website (http://www.ncbi.nlm.nih.gov/geo/). The mentioned databases were searched using the following common terminology: sepsis, human, GPL570 (Affymetrix Human Genome U133 Plus 2.0 Array or hgu133plus2.db), GPL6947 (Illumina HumanHT-12 V3.0 expression beadchip), and GPL10558 (Illumina HumanHT-12 V4.0 expression beadchip) [22, 23]. These platforms were selected because they allow the study of genes related to the immune response. GPL570, GPL6947, and GPL10558 consisted of 58075, 48803, and 47231 probes, respectively. Entering GPL570, GPL6947, and GPL10558 with other terminology, 74 datasets or series were found. Among them, the GSE54514 [11], GSE57065 [12], GSE64456 [13], GSE66099 [14], GSE72829 [15], and GSE95233 [16] datasets were selected because their study designs were suitable to the purpose of this study [19], which included a healthy control group and a sepsis group, with more than 10 samples in each group. Datasets with drug treatment, studied tissue other than whole blood, and known sepsis etiology were excluded (S1 Fig in S1 File). All the datasets used transcriptomes derived from whole blood. Data regarding the time point of gene expression were derived from the first sample after admission. GSE66099 and GSE72829 were composed of several unique datasets, and healthy control and septic shock data from GSE66099 were selected. Healthy control cases and initial sepsis cases were used for each dataset. In the case of GSE72829, a healthy control and a definite bacterial infection group were selected. Randomized selection of sepsis cases without replacement was performed to match the number of healthy control cases and form a balanced correlation matrix. GSE95233 was additionally analyzed for survivor and nonsurvivor sepsis cases (GSE95233-1). Baseline characteristics of the datasets are listed in Table 1.

Table 1

Baseline characteristics of recruited datasets for adult and paediatric patient.

	GSE54514	GSE57065	GSE95233	GSE64456	GSE66099	GSE72829
Patient type	Adult	Adult	Adult	Pediatric	Pediatric	Pediatric
Age, mean (yr) a,b	≥18, 60.3	≥18, 62	≥18, 65	≤2 mo, 29 d	≤10, 3.7	≤17, 1.9 mo
Sex (male %)b	40	67.9	65	52	58	62
Nonsurvivor (%)	25.7	17.9 (28 d)	NA	NA	16.3	19.2
Cohort description	ICU with sepsis or	ICU with	ICU with	ED with sepsis	PICU with sepsis or	Discovery group
	septic shock	septic shock	septic shock	(PECARN study)	septic shock	(IRIS study)
Severity evaluation	APACHEII	SOFA	SOFA	YOSb	PRISM	NA
Control (n)	18	25	22	19	49	52
Sepsis or septic shock (n)	35	28	51	89	198	94
Randomized selection without replacement (septic shock)	18	25	22	-	-	-
Randomized selection without replacement (sepsis)	-	-	-	18	49	52
Measured frequency, frequency (time points)	5 (1, 2, 3, 4, 5 d)	3 (0, 24, 28 hr)	3 (1–2, 3–4, 7–10 d)	1	1	1
Recruited time pointsc	0–24 hr	0–24 hr	1–2 d	NA	0–24 hr	5 d
Platform	Illumina, GPL6947	Affymetrix, GPL570	Affymetrix, GPL570	Illumina, GPL10558	Affymetrix, GPL570	Illumina, GPL10558
Reference	[11]	[12]	[16]	[13]	[14]	[15]

ED, emergency department; ICU, intensive care unit; NA, data not available; hr, hour; d, day; APACHEII, acute physiology and chronic health evaluation II; SOFA, sequential (sepsis related) organ failure assessment.

^aGSE57065 and GSE72829 were summarized using median age.

^b Percent of male and age are summarized data derived from the previous studies, which could be different from the analysis performed in this study. The sepsis patient was randomly selected corresponding to the number of the control cases.

^cYOS, Yale observation score (6–30).

Sepsis diagnosis

Employed gene expression datasets used sepsis diagnosis based on sepsis criteria [1]. The patients were classified as having systemic inflammatory response syndrome if they met two or more of the following criteria: WBC count ≥ 12.0 × 10⁹/L or ≤ 4.0 × 10⁹/L, or presenting more than 10% immature cells; body temperature ≥ 38° C or ≤ 36° C; respiratory rate ≥ 20 breaths per minute, or heart rate ≥ 90 beats per minute. Sepsis was diagnosed with SIRS and documented bacterial infection, and septic shock was diagnosed with sepsis and hypotension not responding to fluid resuscitation [1].

Bioinformatics

For differentially expressed genes (DEG) analysis between control and patient groups based on total transcriptome, gene set enrichment for pathway analysis and network construction, several appropriate R software packages were used [24], including Bioconductor [25], affy [26], GEOquery [27], lumi [28], and limma [29]. The data were downloaded from the GEO database, and the data matrix of gene expression values was formatted using the GEOquery package [27]. In the case of Affymetrix gene expression data sets, the raw microarray data (CEL files) were normalized using the Robust Multichip Averaging (RMA) method, and log2 transformation was performed [25–27]. In the case of Illumina gene expression data sets, the arrays were normalized using the quantile method [27, 29]. The data were annotated with corresponding gene names using R packages hgu133plus2.db, illuminaHumanv3.db, and illumina-Humanv4.db [30, 31].

Gene set enrichment for pathway analysis

Analysis of DEG was performed for total transcriptome. After a t-test using basic functions in R software, genes were selected based on false discovery rate of less than 0.01 (FDR < 0.01) [32]. Based on DEG analysis, gene set enrichment for pathway analysis using GAGE package, based on R language, was performed [33]. GAGE package applies input data to the Kyoto Encylopedia of Genes and Genomes (KEGG) pathways. Significant pathways were plotted for heatmap that was selected based on global P-value of less than 0.05.

Immune response related genes

To analyze the function of the immune response in the control and sepsis groups, KEGG immune system pathways were used [34]. Twenty-one pathways of the immune system were included, as follows: hematopoietic cell lineage (has04640); complement and coagulation cascade (has04610); platelet activation (hsa04611); Toll-like receptor signaling (hsa04620); Toll and Imd signaling (hsa04624); NOD-like receptor signaling pathway (hsa04621); RIG-I-like receptor signaling (hsa04622); cytosolic DNA-sensing (hsa04623); C-type lectin receptor signaling (hsa04625); natural killer cell-mediated cytotoxicity (hsa04650); antigen processing and presentation (hsa04612); T cell receptor signaling (hsa04660); Th1 and Th2 cell differentiation (04658); Th17 cell differentiation (hsa04659); B cell receptor signaling (hsa04662); Fc epsilon RI signaling (hsa04664); Fc gamma R-mediated phagocytosis (hsa04666); leukocyte transendothelial migration (hsa04670); intestinal immune network for IgA production (hsa04672), and chemokine signaling (hsa04062). The 21 pathways were assigned to adaptive immunity, innate immunity, NK cell activity, leukocyte migration, antigen presentation, complement-coagulation, platelet activity, cytokines-chemokines, signaling, and hematopoiesis. From them, 1908 genes in these pathways and 998 unique genes were extracted (S1 and S2 Tables in S1 File).

As genes were included in multiple pathways, we utilized Gene Ontology (GO) and Reactome pathways for additional information [35, 36]. From the GO pathway, we included activation of the innate immune response, complement activation, and immune response-activating signal transduction. From the immune response pathway, the adaptive immune response and the innate immune response were included. Leukocyte homeostasis, leukocyte activation, and leukocyte migration were also collected. From the Reactome database, we included the hemostasis pathway; platelet homeostasis; platelet adhesion to exposed collagen; platelet activation, signaling, and aggregation; formation of fibrin clots; dissolution of fibrin clots, and cell surface interactions at the vascular wall. In the immune system pathway, adaptive immune system, innate immune system, and cytokine signaling in the immune system were included (S3 Table in S1 File).

Genes included in multiple pathways were assigned to one pathway with their most frequent function among the KEGG, GO, and Reactome pathways. If the frequency of the gene that was included in the pathway was equal, undetermined (U) value was assigned. If the frequency was equal and included chemokines-cytokines, signaling, or hematopoiesis pathways, pathways other than these ones were assigned (S3 Table in S1 File).

Network topology analysis

The network was constructed using Spearman’s correlations equal to or greater than 0.8 and a P-value of less than 0.05 for each set of molecules. Each gene was regarded as a node, and the correlation pairs between the genes were regarded as links [37]. Links were created between the nodes with significant correlation coefficients. Network analysis results were plotted, and the topologic parameters were calculated using Cytoscape and an additional application, NetworkAnalyzer [37–39]. The clustering coefficient is defined as the number of links between associated nodes. The degree of a node is defined as the number of links connected to a node (gene), and a hub is defined as a node with a relatively high degree. Network density is the ratio of links in the network to the maximum possible number of links. The length of a path is the number of unique links between the nodes that are present in the path. Distance is the shortest path length between two nodes. The network diameter is the maximum length of the shortest path. The characteristic path length is the average shortest path length and provides the expected distance between two nodes. The average number of neighbors is the number of identical nodes connected to a node. The normalized version of the average number of neighbors is the network density. Network heterogeneity is the variance in the number of links of a node divided by the number of mean links. The average number of neighbors indicates the average connectivity of the nodes. The topological coefficient is a relative measure of the extent to which a node shares neighbors with other nodes. Network centralization is close to one when the network resembles a star-like figure, while uniformly connected networks have a centralization value close to zero. The stress centrality of a node is the number of shortest paths passing through that node. The betweenness centrality of a node is the number of shortest paths between two other nodes that pass through that node, divided by the total number of shortest paths between the two nodes; betweenness centrality reflects the amount of control that this node exerts over other nodes. Closeness centrality reflects how close nodes are, which implies that nodes with high closeness centrality values have the shortest distance to other nodes on average. Eccentricity reflects the maximum distance of a node from other nodes [38, 39]; a high eccentricity value means that all other nodes are in proximity, a low eccentricity value means that there are at least one node and its neighbors that are far from that node. Degree correlation defines the relationship between the degree of a node and that of connected nodes [40]. Modularity is defined as the tendency of nodes to cluster and to interact, coordinating biological processes [41, 42]. The parameters were analyzed using comparison of proportions (MedCalc Software version 18.2.1, Mariakerke, Belgium).

Gene set enrichment for pathway analysis

Total transcriptome analysis of healthy control and sepsis group revealed that 176 to 8229 genes showed expression differences in datasets. Gene set enrichment for pathway analysis was performed based on KEGG pathways (S4 Table in S1 File). Among up-regulated or down-regulated pathways with statistical significance, we searched for pathways associated with immune-related pathway (S2–S6 Figs in S1 File). We found that complement and coagulation cascade pathway and Fc epsilon RI signaling and Fc gamma R-mediated phagocytosis pathways were up-regulated in 2/3 of datasets from adult patients and in 3/3 of those from pediatric patients. T cell receptor signaling pathway was down-regulated in 2/3 and 3/3 of datasets from adult and pediatric patients, respectively. Antigen presentation and processing were down-regulated in 2/3 of datasets from adult and pediatric patients (Table 2). Altogether, these findings were interpreted along with the network analysis results. In network analysis, adaptive immunity genes were prominent and isolated from network in sepsis patient group (Fig 1). T cell receptor signaling (TCR) pathway was down-regulated in 5/6 of the datasets. These data imply that adaptive immunity related with T cells might be impaired or dysregulated in sepsis. In addition, adaptive immunity was associated with unfavorable prognosis in sepsis group. Complement and coagulation cascade was prominent in sepsis survivors of GSE95233 datasets (Fig 2). Complement and coagulation cascade was increased in sepsis in 5/6 of datasets and these data imply that these pathways may be involved in convalescence process.

Fig 1

Network structure of GSE57065 of the (A) control group and (B) sepsis group. In the sepsis group, the adaptive immunity network was prominent and isolated (green).

Fig 2

Network structure of GSE95233 of (A) sepsis survivors and (B) sepsis nonsurvivors. In the survivor group, the platelet (olive green) and coagulation cascade network (yellow) were isolated and prominent. In the nonsurvivor group, the adaptive immunity (green) and innate immunity networks (blue) were prominent and isolated.

Table 2

Differentially expressed gene analysis between control and sepsis followed by gene set pathway analysis^a.

	GSE54514	GSE57065	GSE95233	GSE64456	GSE66099	GSE72829
Patient type	Adult	Adult	Adult	Pediatric	Pediatric	Pediatric
Control cases (n)	18	25	22	19	49	52
Sequential selection of sepsis cases (n)	18	25	22	19	49	52
Differentially expressed Genes (n)	176	7885	6784	3123	8229	7463
Up regulated pathway (n)	0	60	44	13	78	51
Down regulated pathway (n)	0	33	33	16	33	15
Up regulated immune pathway		C, Fc, RIG	C, Fc	C, Fc, NK, K, TLR	APC, C, CXC & etc, Fc, IgA,	BCR, C, Fc, CXC & etc, LK, NK, TLR
Down regulated immune pathway		APC, TCR, H, K	APC, BCR, CXC & etc, IgA, TCR	APC, IgA, TCR	H, TCR	APC, IgA, TCR

^aDifferentially expressed gene was selected based on FDR (P<0.01).

APC, Antigen processing and presentation; BCR, B cell receptor signaling; C, Complement and coagulation; CXC & etc, chemokine signaling; Fc, Fc gamma R mediated phagocytosis or Fc epsilon RI signaling; H, hematopoietic cell lineage; IgA, intestinal immune network for IgA production; LM, leukocyte endothelial migration; NK, NK cell mediated cytotoxicity; RIG, RIG-I like receptor signaling; TLR, Toll like receptor signaling; TCR, T cell receptor signaling.

Network structure and immune function

The network topologic structure revealed that, in sepsis, adaptive immunity genes (green colored) formed a prominent cluster, isolated from the rest of the network (Fig 1 and S7–S18 Figs in S1 File). In that isolated cluster, CD247 (CD3 Zeta chain), CD8A, ITK (tyrosine protein kinase ITK / TSK), LAT (linker for activation of T cells), LCK (leukocyte C-terminal Src kinase) were found in all 6 datasets. CD2, FYN (Src family tyrosine kinase), GATA3 (GATA-binding protein 3), IL7R, RASGRP1 (RAS guanyl-releasing protein 1) were all found in adult datasets, but 2/3were found in the pediatric group. CBLB (E2 ubiquitin protein ligase CBL-B), CD3D, CD3G, ZAP70 (Zeta-chain-associated protein kinase 70) were found in 2/3 of adult datasets, but all were found in pediatric datasets. Other genes are listed in S5 Table in S1 File. In addition, NK cell activity (light blue) tended to be included in the network in the sepsis group, whereas it was isolated in the healthy control group. The nonsurvivor network from GSE95233 revealed that adaptive immunity was isolated and prominent from the rest of the network, whereas platelet activity (olive green) and the compliment and coagulation cascade pathway (yellow) were isolated and prominent in the survivor group (Fig 2). Platelet, complement and coagulation-related genes were prominent in survivors; C6, C7, F11, COL1A1, COL3A1, ADCY6, GP5, FGG, MYLK4, SERPINE1, SERPIND1, CR2, etc. were included. CD8B, IL1A, IL3, IL12B, IL21, IL23R, and other genes were also included. The isolated gene cluster from nonsurvivors included CD2, CD3D, CD3G, CD8A, CD247, CBLB, GATA3, ITK, LCK, CTLA4, IL10, ZAP70, and other genes.

Network parameter analysis

Comparison of topologic parameters revealed that the clustering coefficient in sepsis was increased in 2/3 of adult datasets and in 3/3 of pediatric datasets (Table 3). Network heterogeneity was increased in 3/3 of adult datasets, but was decreased in 3/3 of pediatric datasets. Modularity value was increased in 1/3 of adult datasets, whereas modularity value was increased in 3/3 of pediatric datasets in the sepsis group. Shortest path was decreased in 2/3 of adult and pediatric datasets, whereas average degree was increased in 2/3 of adult and pediatric datasets (Table 3). For GSE95233-1 analysis, difference between network parameters was calculated as follows: (case—control / control) x 100. The difference of clustering coefficient between survivor (0.390) and nonsurvivor (0.337) was -13.5%. The difference of modularity between survivor (0.581) and nonsurvivor (0.694) was 19.4%.

Table 3

Network topological parameters from gene expression data from GEO dataset for adult and paediatric patient.

	GSE54514		GSE57065		GSE95233		GSE95233-1		GSE64456		GSE66099		GSE72829
Patient type	Adult		Adult		Adult		Adult		Pediatric		Pediatric		Pediatric
Disease type	Control	septic shock	control	septic shock	control	septic shock	survivor	non survivor	Control	sepsis	Control	sepsis	Control	sepsis
Clustering coefficient	0.319	0.323	0.328	0.280	0.328	0.367	0.390	0.337	0.367	0.390	0.337	0.388	0.381	0.435
Network density	0.013	0.016	0.012	0.010	0.012	0.039	0.026	0.014	0.039	0.026	0.014	0.012	0.030	0.034
Network heterogeneity	1.039	1.142	1.151	1.172	1.151	1.504	1.358	1.193	1.504	1.358	1.193	1.091	1.249	1.129
Connected components	27	30	42	61	42	34	33	47	34	33	47	32	16	25
Network diameter	19	17	14	15	14	12	19	11	12	19	11	14	13	9
Network centralization	0.059	0.079	0.103	0.048	0.103	0.196	0.113	0.075	0.196	0.113	0.075	0.078	0.148	0.121
Shortest path	160514	80598	59638	17222	59638	142640	121152	11730	142640	121152	11730	36902	38092	22532
Characteristic path length	7.418	5.267	5.276	5.088	5.276	3.862	7.211	3.482	3.862	7.211	3.482	4.705	3.92	2.824
Average degree, <k>	5.93	7.42	5.30	3.64	5.30	17.87	11.71	4.04	17.87	11.71	4.04	4.50	8.43	10.27
Number of nodes, N	464	453	426	367	426	457	448	297	457	448	297	380	284	299
Degree correlation, μ	0.489	0.508	0.422	0.810	0.422	0.403	0.715	0.546	0.403	0.715	0.546	0.210	0.475	0.346
Modularity, Mc	0.743	0.662	0.644	0.806	0.644	0.198	0.581	0.694	0.198	0.581	0.694	0.745	0.434	0.603

GSE95233-1 indicates survivors (n = 16) and non-survivors (n = 16) among sepsis patients.

We searched for nodes with such a high degree that they could be considered as hub nodes. In the GSE54514 control and sepsis groups, CD3D and ITPR3 (inositol 1,4,5 trisphosphate receptor type 3) had degree 33 and 43, respectively. In the GSE57065 control and sepsis groups, ATG5 (autophagy-related 5), and GBP5 (guanylate-binding protein 5) had degree 19 and 25, respectively. In the GSE95233 control and sepsis groups, IFNGR1 (interferon gamma receptor 1), and ESAM (endothelial cell selective adhesion molecule) had degree 49 and 21, respectively.

In the GSE64456 control and sepsis groups, LYN (tyrosine protein kinase Lyn), and IFI16 (interferon gamma-inducible protein 16) had degree 56 and 70, respectively. In the GSE66099 control and sepsis groups, NCF2 (neutrophil cytosolic factor 2), NLRP12 (NLR family pyrin domain-containing 12), and NCR2 (natural cytotoxicity-triggering receptor 2) had degree 26, 26, and 34, respectively. In the GSE72829 control and sepsis groups, C5AR1 (complement C5a receptor 1) and CSF3R (colony-stimulating factor 3 receptor) had degree 50 and 46, respectively. These network parameters are listed in S6–S17 Tables in S1 File.