This article is mentioned in:

Introduction

Sepsis is defined as the systemic inflammatory response to infection and is one of the leading causes of mortality in critically ill patients despite the application of numerous antibiotics and resuscitation therapies (1). Sepsis syndromes could be classified as systemic inflammatory response syndromes (SIRSs), sepsis, severe sepsis and septic shock according to the American College of Chest Physicians and the Society of Critical Care first published Care (2). In addition, the incidence of sepsis is increasing worldwide. In the USA the current incidence of sepsis is ~3 in 1,000 people, whereas severe sepsis leads to at least 200,000 deaths per year (3). Furthermore, severe sepsis and septic shock account for 30–50% of hospital-reported mortality (4). Recently, consensus on the treatment of sepsis includes advanced supportive care in the intensive care unit and use of bundle therapies (5). However, due to the non-specific nature of the signs and symptoms of sepsis, the diagnosis and treatment of sepsis are complicated.

Recently, numerous studies have been performed to identify the pathogenesis of sepsis (6–8). A number of biomarkers can be used in the diagnosis of sepsis, however, none of them has sufficient specificity or sensitivity in the clinical setting (9–11). C-reactive protein (CRP) and procalcitonin (PCT) have been widely used because of their relatively better specificity and prognostic capability (12–14). The concentration of CRP is <0.8 mg/l and can increase 1,000-fold in response to an acute-phase stimulus (15). CRP is a protein that is synthesized in the liver and rises in response to inflammation (16). Moreover, it may help macrophages remove microorganisms by binding the phospholipid components (17). Nowadays, CRP is also treated as a biomarker for evaluating sepsis severity and prognosis or to monitor treatment response (12). PCT, which is produced by parafollicular cells of the thyroid and neuroendocrine cells of the lungs and the intestines is a 116 amino acid polypeptide precursor of the hormone calcitonin. PCT was first linked to infectious disease by Assicot et al (18) and was formally proposed as an adjunctive diagnostic biomarker in 2008 (19). It is maintained at a low level in healthy people and increases 1,000-fold during active infection (20). Furthermore, there are also several meta-analyses demonstrating that PCT could be used as a diagnostic marker in sepsis (21,22).

Nevertheless, there is currently no gold biomarker that exists as a marker of sepsis. Thus, identification of a new biomarker is urgently required. In order to further identify the molecular pathogenesis of sepsis, microarray data were firstly downloaded, then the raw data were analyzed to construct a protein-protein interaction (PPI) network. Subsequently, differentially expressed clusters in the PPI network were identified and significantly enriched pathways and functions of the genes in the clusters were also screened. Finally, potential molecular markers were identified using the support vector machine (SVM) method.

Materials and methods

Obtaining and preprocessing of mRNA expression profile data

The mRNA expression profiles of sepsis and non-sepsis samples were obtained from the National Center of Biotechnology Information Gene Expression Omnibus database. The access number was GSE12624 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE12624) and the datasets of 36 samples with septic shock following trauma (sepsis samples) and 34 samples without septic shock following trauma (non-sepsis samples) were used for further analysis. The platform used here was called GPL4204 GE Healthcare/Amersham Biosciences CodeLink UniSet Human I Bioarray. The original data at probe symbol level were first converted into expression values at gene symbol level. Next, missing data was imputed and median data normalization was performed using robust multichip averaging (23). Besides, principal component analysis (PCA) (24), which was used as a computational procedure for biomarker identification and for the classification of multiclass gene expression was performed to identify the difference between sepsis and non-sepsis samples.

PPI network construction

PPIs illustrate valuable information for the elucidation of cellular function, and protein interaction studies have been developed to be a focal point of recent biomolecular research. The Human Protein Reference Database (HPRD) (25) is a novel protein information resource illustrating various features of proteins, including the domain architecture, molecular function, tissue expression, subcellular localization, enzyme-substrate correlation and PPIs. In the present study, all the human PPI pairs in HPRD were initially collected. Next, the Pearson correlation coefficients for all the interacting genes were calculated based on their expression values under the sepsis and non-sepsis status with a coefficient <0.5 used as the cut-off criterion. This was done to obtain the PPI networks under these two statuses. Furthermore, Cytoscape (26) was used to visualize the PPI networks in order to further observe the correlation between genes.

Hierarchical clustering and analysis of covariance (ANCOVA) global test for differentially expressed clusters

Hierarchical clustering is a method of cluster analysis that seeks to build a hierarchy of clusters (27). Euclidean distance was selected as a measure of distance between pairs of genes in the PPI network. The present study used the package hclust (http://CRAN.R-project.org/package=gplots) in R language to perform the hierarchical clustering of two PPI networks, with the requirement that each cluster should have had >5 genes. Finally, Package GlobalAncova in R language was employed to identify the differentially expressed cluster using the ANCOVA global test (28), which focuses on phenotype effects and gene-phenotype interactions. P<0.05 was defined as a threshold.

Enrichment analysis of differentially expressed clusters

In order to study differentially expressed clusters at the functional level, Gene Ontology (GO) functional enrichment (29) and Kyoto Encyclopedia of Gene and Genomes (30) pathway enrichment analyses were performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) software (31). DAVID software has been widely used to identify biological processes involving a given list of genes. In the present study, a fold change discovery of ≤0.05 was set as the cut-off criterion for enrichment analysis.

Identification of molecular markers by the SVM method

The SVM method has been demonstrated to be an useful classification and regression method that uses machine learning theory to maximize the predictive accuracy while avoiding overfitting of data (32). Differential clusters were initially ranked according to their P- and F-values (F test), and the cluster with the highest P-value was then defined as the class feature. Secondly, the SVM method was employed to classify the samples using the ksvm function in the kernlab package in R language, and the 10-fold cross-validation method was performed to evaluate the classification results. A feature selection error rate of <0.1 was selected as the criteria.

Results

Preprocessing of mRNA expression profile data

A total of 7,672 genes were obtained after preprocessing of the mRNA expression profile data. PCA analysis based on gene expression values revealed that sepsis samples had a great similarity with non-sepsis samples (Fig. 1). Therefore, further bioinformatics analyses were conducted to identify the molecular markers for distinguishing sepsis samples from non-sepsis samples.

Figure 1. Principal component analysis of gene expression in sepsis and non-sepsis samples. Blue points represent the sepsis and red points represent non-sepsis samples. PC1, principal component 1.

PPI network construction

The HPRD database was used to construct the PPI network and Cytoscape software was employed to visualize the network. The PPI network of genes under sepsis collected 1,996 genes and 2,645 interactions between them; the PPI network of genes under non-sepsis status collected 2,147 genes and 2,783 interactions (Fig. 2). Further analysis revealed that there were 1,438 overlapping genes and only 992 overlapping interactions between sepsis and non-sepsis samples (Fig. 3).

Figure 2. Protein-protein interaction network. (A) Sepsis and (B) non-sepsis samples.

Figure 3. Overlapping conditions of two protein-protein networks. (A) Sepsis and (B) non-sepsis samples.

Screening of differentially expressed clusters

In total, 40 clusters with 20 clusters in both PPI networks were identified by hierarchical clustering (Fig. 4). Furthermore, 24 differentially expressed clusters, each with 12 clusters in both networks, were identified by the ANCOVA global test with a threshold of P<0.05.

Figure 4. Hierarchical clustering analysis of clusters in two protein-protein networks. (A) Non-sepsis and (B) sepsis samples.

Enrichment analysis of differentially expressed clusters

Function enrichment analysis revealed that the genes in the 24 differentially expressed clusters were mainly enriched in the following GO terms: Positive regulation of macromolecule metabolic processes, transcription factor binding, regulation of cell death and the intracellular signaling cascade (Table I). However, no pathways were enriched.

Table I. GO analysis of differentially expressed clusters (top 15).

Table I.

GO analysis of differentially expressed clusters (top 15).

Category	Term	FDR
GOTERM_BP_FAT	GO:0010604~positive regulation of macromolecule metabolic process	6.58×10−29
GOTERM_MF_FAT	GO:0008134~transcription factor binding	1.09×10−27
GOTERM_BP_FAT	GO:0010605~negative regulation of macromolecule metabolic process	2.52×10−23
GOTERM_BP_FAT	GO:0010941~regulation of cell death	1.83×10−22
GOTERM_BP_FAT	GO:0007242~intracellular signaling cascade	2.59×10−22
GOTERM_BP_FAT	GO:0043067~regulation of programmed cell death	8.30×10−22
GOTERM_MF_FAT	GO:0030528~transcription regulator activity	7.43×10−21
GOTERM_BP_FAT	GO:0042981~regulation of apoptosis	1.21×10−20
GOTERM_MF_FAT	GO:0004672~protein kinase activity	1.08×10−19
GOTERM_BP_FAT	GO:0016310~phosphorylation	7.73×10−19
GOTERM_BP_FAT	GO:0006468~protein amino acid phosphorylation	8.91×10−19
GOTERM_BP_FAT	GO:0044093~positive regulation of molecular function	4.92×10−17
GOTERM_BP_FAT	GO:0010628~positive regulation of gene expression	5.97×10−17
GOTERM_BP_FAT	GO:0006796~phosphate metabolic process	6.14×10−17
GOTERM_BP_FAT	GO:0006793~phosphorus metabolic process	6.14×10−17

[i] GO, gene ontology; BP, biological process; FAT, functional annotation tool; MF, molecular function; FDR, false discovery rate.

Identification of molecular markers by the SVM method

The 24 clusters that were sorted by P-value are listed in Table II. The SVM method was employed in order to classify the samples according to the cluster sequences listed in Table II. Notably, when the seventh cluster was added, the feature selection error rate of 10-fold cross-validation was >0.1. Therefore, 207 genes in the top six clusters were selected as potential biomarkers, including CDC42, CSF3R, GCA, HMGB2, RHOG, SERPINB1 and TYROBP in cluster 15, as well as SERPINA1, FCER1G and S100P. Fig. 5 illustrates the SVM calculating process. In order to validate the non-randomness of these 207 genes, 207 genes were randomly selected from the mRNA expression profile data 1,000 times and the results of the 10-fold cross-validation were always >0.1. This proved the effectiveness of the selected molecular markers. Further enrichment analysis demonstrated that these 207 genes were mainly enriched in functions of the intracellular signaling cascade, regulation of programmed cell death, regulation of cell death, regulation of apoptosis and leukocyte activation (Table III); whereas no pathways were enriched.

Figure 5. Trend graph of error rate for feature selection using support vector machine analysis method.

Table II. Sorted differentially expressed clusters of sepsis and non-sepsis samples.

Table II.

Sorted differentially expressed clusters of sepsis and non-sepsis samples.

Clusters	F-value	p.perm	p.approx	State
15	9.12	<0.01	2.74×10−6	s
18	6.66	<0.01	4.87×10−3	ns
14	6.03	<0.01	1.06×10−4	s
10	5.26	<0.01	7.52×10−4	ns
13	4.83	<0.01	1.25×10−5	ns
8	4.61	<0.01	1.40×10−3	s
13	3.94	<0.01	4.86×10−6	s
3	3.77	<0.01	2.92×10−6	ns
2	3.58	<0.01	6.30×10−4	s
11	3.25	<0.01	2.37×10−3	ns
3	3.23	<0.01	2.26×10−3	s
11	3.03	<0.01	3.31×10−3	s
2	2.94	<0.01	6.03×10−3	ns
5	2.82	<0.01	2.06×10−2	s
6	2.82	<0.01	6.06×10−3	s
4	2.64	<0.01	8.88×10−3	ns
16	4.88	0.01	1.11×10−3	ns
6	3.76	0.01	6.36×10−4	ns
4	3.42	0.01	2.21×10−3	s
8	2.94	0.01	5.65×10−3	ns
9	2.43	0.02	4.01×10−2	s
12	2.26	0.02	3.19×10−2	s
1	2.65	0.03	1.53×10−2	ns
5	2.36	0.04	1.87×10−2	ns

[i] s, septic samples; ns, non-septic samples; p-perm, p-values from the permutation test; p.approx, P-values by means of an approximation for a mixture of χ² distribution.

Table III. GO analysis of molecular markers (top 15).

Table III.

GO analysis of molecular markers (top 15).

Category	Term	FDR
GOTERM_CC_FAT	GO:0005829~cytosol	1.33×10−8
GOTERM_BP_FAT	GO:0007242~intracellular signaling cascade	3.76×10−5
GOTERM_CC_FAT	GO:0031982~vesicle	9.54×10−3
GOTERM_CC_FAT	GO:0044459~plasma membrane part	1.12×10−2
GOTERM_BP_FAT	GO:0010033~response to organic substance	1.39×10−2
GOTERM_BP_FAT	GO:0043067~regulation of programmed cell death	1.40×10−2
GOTERM_CC_FAT	GO:0005886~plasma membrane	1.46×10−2
GOTERM_BP_FAT	GO:0010941~regulation of cell death	1.50×10−2
GOTERM_CC_FAT	GO:0015629~actin cytoskeleton	1.67×10−2
GOTERM_MF_FAT	GO:0032403~protein complex binding	2.73×10−2
GOTERM_BP_FAT	GO:0042981~regulation of apoptosis	3.31×10−2
GOTERM_BP_FAT	GO:0045321~leukocyte activation	6.04×10−2
GOTERM_BP_FAT	GO:0031400~negative regulation of protein modification process	7.36×10−2
GOTERM_CC_FAT	GO:0009986~cell surface	7.38×10−2
GOTERM_CC_FAT	GO:0031988~membrane-bounded vesicle	7.53×10−2
GOTERM_BP_FAT	GO:0016192~vesicle-mediated transport	9.46×10−2

[i] GO, gene ontology; CC, cellular components; FAT, functional annotation tool; BP, biological process; MF, molecular function; FDR, false discovery rate.

Discussion

Sepsis and its complications are a common cause of infectious disease and hospital-reported mortality worldwide (33). The present study aimed to investigate the potential mechanism of sepsis, and to identify genes that can be used for diagnosing and developing candidate molecularly targeted therapy. In total, 7,672 genes were obtained after preprocessing of the mRNA expression profile data. Following hierarchical clustering analysis and the ANCOVA global test, 24 differentially expressed clusters with 12 clusters in each PPI network were identified. Moreover, 207 genes in the top six clusters were selected using SVM, and the functional enrichment analysis revealed that they were mainly enriched in the intracellular signaling cascade, regulation of programmed cell death, regulation of apoptosis and leukocyte activation.

The ANCOVA global test identified 24 differentially expressed clusters, and cluster 15 had the highest P-value. Genes in this cluster, including CDC42 (34), CSF3R (35), GCA (36), HMGB2 (37), RHOG (38), SERPINB1 (39) and TYROBP (40) had already been linked with sepsis. SERPINA1, FCER1G and S100P were also genes in this cluster. Gene SERPINA1 encodes alpha-1-antitrypsin, which is a serine protease inhibitor. Moreover, the targets of SERPINA1 include elastase, plasmin, thrombin, trypsin, chymotrypsin and plasminogen activator, which participates in inflammatory processes (41). Buttenschoen et al (42) demonstrated that the diagnostic value of SERPINA1 levels could be applied in order to distinguish sepsis from SIRS and to assess prognosis. Recently, Su et al (43) revealed that SERPINA1 was downregulated in patients with sepsis compared with SIRS patients, and further analysis demonstrated that SERPINA1 was involved in sepsis differentiation. The FCER1G gene encodes the γ-subunit of Fc epsilon RI (FcRγ), which is an immunoreceptor tyrosine-based activation motif-bearing signal transduction subunit of the Fc receptor family (44). The FCER1G gene was upregulated in sepsis according to Hu et al (45). Furthermore, it has a deleterious effect on sepsis, and FcRγ-/− mice demonstrated an increased survival during sepsis due to increased Escherichia coli phagocytosis (46). The S100P gene encoding the S100 calcium binding protein, which is a member of the S100 family proteins, contains 2 EF-hand calcium-binding motifs. Sepsis-associated encephalopathy (SAE) is the organ dysfunction accompanied with sepsis (47). S100P is produced mainly by the central nervous system, and the elevated serum level of S100P is a biomarker of neuronal damage occurring in SAE (48). Therefore, elevated serum levels of S100P may be associated with critical illness and may be treated as the biomarkers of brain damage during sepsis.

GO functional enrichment analysis revealed that 207 genes identified by SVM, including TLR2 and RAB27A, were mainly enriched in the intracellular signaling cascade, regulation of programmed cell death and cell death, regulation of apoptosis and leukocyte activation. Furthermore, programmed cell death is an important mechanism during the immunopathogenesis of sepsis. Notably, apoptosis is one form of programmed cell death. In addition, early programmed cell death of lymphocytes destroys innate and adaptive immunity, which would reduce the ability of protecting against pathogens (49). Also, programmed cell death of parenchymal cells in the lung, liver and gut would facilitate organ failure and death (50). The TLR2 gene encodes Toll-Like Receptor 2, which is a member of the Toll-like receptor family expressed on the macrophage recognizing pathogen-associated molecular patterns (51). Several reports have suggested the role of TLR2 in the induction of pathogen-induced programmed cell death (52–54). The present study identified that TLR2 demonstrated differential expression in sepsis samples, which was consistent with the observations of Armstrong et al (55). Therefore, we inferred that TLR2 may be involved in sepsis by interrupting programmed cell death. The RAB27A gene encoding guanosine triphosphate (GTP)-binding protein Ram belongs to the GTPase superfamily, Rab family (56). Ménasché et al (57) have reported that genetic defects in Rab27a may lead to immunodeficiency in humans caused by programmed cell death. Furthermore, Johnson et al (58) illustrated that Rab27a deficiency is associated with increased survival and reduced neutrophil infiltration of the liver in a model of lipopolysaccharide-induced systemic inflammation. Therefore, Rab27a may participate in organ failure which accompanies sepsis.

Overall, in order to illustrate the pathological mechanisms underlying sepsis, gene expression profiles containing 70 samples were downloaded and analyzed. SERPINA1, FCER1G and S100P in the selected differential clusters may be potential biomarkers. Moreover, TLR2 and Rab27a may exert certain roles in sepsis by interrupting programmed cell death. However, more experiments are required in order to confirm these results.

References