Co-expressed gene modules identified in sepsis
Although sepsis transcriptome has been profiled in several studies, the small sample size and microarray platform limit the analysis power. We used the currently largest sepsis RNA-Seq dataset GSE185263 to construct a gene co-expression network (Fig. 1A). Using the WGCNA method, a total of 37 co-expressed gene modules were identified (Table 2). To test if these modules are stable, the module stability analysis was performed. Results showed that the average module stability was larger than 0.96 with a standard deviation of 0.015 (Fig. 1B). Functional enrichment analysis revealed that these modules were associated with transcription factors, translation, immune response, cell cycle, mitochondrion, leukocyte activation, platelet activation, and interferon signaling (Table 2). All the module genes and related information is provided in Supplementary Table 1.
Table 2 Functional and immune cell annotation of the identified modules
Module (size)
|
Functional annotation
|
Hub gene
|
Immune cell
|
1 (1784)
|
Factor: hdac2 (1E-56) hsa-miR-21-5p (1E-35) nucleus (1E-29)
|
GABPA
|
Th1 CD4+ T cell (-)
|
2 (2572)
|
Factor: hdac2 (7E-41) Nucleic acid metabolism (1E-37) nucleus (2E-33)
|
PDS5A
|
CD4+ T cell
|
3 (1592)
|
Factor: SP2 (6.4E-58) Factor: ZF5 (1E-57)
|
RABL6
|
Hematopoietic stem cell (-)
|
5 (1438)
|
Factor: ETF (5E-18) hsa-miR-21-5p (3E-10) whole membrane (3E-9)
|
SAMD8
|
Neutrophil;
Th1 CD4+ T cell (-)
|
6 (919)
|
Factor: Sp1 (1E-34) Factor: ETF (1E-32)
|
MAPK3
|
NKT; Memory CD4+ T cell (-)
|
14 (993)
|
Generation of second messenger molecules (1E-20) Adaptive immune response (2E-16)
|
DDX24
|
CD8+ T cell
|
17 (248)
|
Cell cycle (1E-87) hsa-miR-193b-3p (3E-37) Kinetochore (1E-36)
|
NUSAP1
|
Th2 CD4+ T cell; Neutrophil (-)
|
18 (245)
|
Factor: Elf-1 (0.006)
|
RNU7-181P
|
|
20 (1160)
|
Factor: ER81 (7E-38) Factor: ZF5 (1E-36)
|
COPE
|
Hematopoietic stem cell (-)
|
23 (377)
|
Erythrocyte differentiation (7E-7)
|
TRIM10
|
|
27 (160)
|
antigen binding (1E-171) Adaptive immune response (6E-121)
|
IGHV3-21
|
Plasma B cell
|
30 (290)
|
Mitochondrion (5E-28) Mitochondrial translation (2E-20)
|
SNRPC
|
Th1 CD4+ T cell; Neutrophil (-)
|
33 (1139)
|
Myeloid leukocyte activation (2E-18)
|
C3orf86
|
Neutrophils; CD8+ T cell (-)
|
34 (225)
|
Defense response to virus (9E-40) Interferon Signaling (7E-34) Innate immune response (5E-33)
|
RTP4
|
Neutrophil
|
35 (151)
|
Mitochondrial inner membrane (1E-26) Factor: ER71 (2E-15)
|
ERH
|
Common lymphoid progenitor
|
36 (112)
|
Factor: SP1 (1E-7)
|
UBL7
|
|
37 (166)
|
Platelet activation, signaling and aggregation (6E-18)
|
ABLIM3
|
|
39 (101)
|
Factor: Erg (2E-6) Organelle lumen (9E-6)
|
ZC3H15
|
Common lymphoid progenitor
|
40 (169)
|
ncRNA metabolism (5E-15)
|
DDX1
|
Common lymphoid progenitor; Neutrophil (-)
|
42 (92)
|
Formation of a pool of free 40S subunits (2E-91) Protein targeting to ER (4E-86) Viral mRNA Translation (4E-84)
|
RPS3
|
Neutrophil (-)
|
43 (89)
|
B cell receptor signaling pathway (1E-5)
|
PAX5
|
B cell; Monocyte (-)
|
46 (82)
|
Plasma membrane (0.002)
|
ZNF385A
|
|
49 (65)
|
Exocytosis (2E-19) Neutrophil activation involved in immune response (5E-19)
|
AQP9
|
Neutrophil; B cell (-)
|
50 (60)
|
Chromatin organization (3E-7) hsa-miR-218-5p (4E-6)
|
ARID1A
|
|
52 (58)
|
Natural killer cell mediated cytotoxicity (3E-6) Factor: MAFA (0.02)
|
TBX21
|
|
55 (58)
|
hsa-miR-484 (3E-5) Factor: CTC (8E-5)
|
MYO9B
|
Common lymphoid progenitor
|
56 (316)
|
Neutrophil activation (2E-12)
|
RASGRP4
|
Neutrophil; Central memory CD8+ T cell (-)
|
57 (51)
|
Chromatin organization (5E-12)
|
YLPM1
|
CD4+ T cell
|
59 (49)
|
Neutrophil degranulation (3E-23) Neutrophil activation (5E-22)
|
CEACAM6
|
Common myeloid progenitor; CD8+ T cell (-)
|
60 (47)
|
|
|
|
61 (87)
|
Factor: ER81 (7E-6)
|
LAMTOR3
|
|
63 (45)
|
|
SLC51A
|
Endothelial; CD4+ T cell (-)
|
64 (45)
|
Peptide disulfide oxidoreductase activity (0.006)
|
NDUFB3
|
CD4+ T cell (-)
|
66 (38)
|
Eukaryotic translation elongation (1E-37) Selenocysteine synthesis (7E-35) Influenza Infection (4E-30)
|
RPL31
|
|
67 (36)
|
hsa-miR-218-2-3p (0.01)
|
USP12
|
|
69 (34)
|
Interferon alpha/beta signaling (6E-7)
|
RBCK1
|
pDCs
|
70 (31)
|
Ribosome biogenesis in eukaryotes (0.01)
|
LEO1
|
Neutrophil (-)
|
pDCs: Plasmacytoid dendritic cells (pDCs); (-): the immune cell proportion is negatively associated with module expression.
Modules were associated with immune cell populations
To check if these identified modules were associated with immune cell proportion, CIBERSORT was used to deconvolute GSE185263 into immune cell proportion. Correlation analysis between immune cell proportion and module expression indicates the links between them (top assignment with an absolute R >0.6, Table 2).
Modules were differentially expressed in sepsis and control
As the high dimensional sepsis transcriptome data has been reduced to tens of modules, module expression was compared between different sepsis groups. We found 27 modules were differentially expressed between sepsis and normal groups. Among them, 17 modules were up-regulated and 10 modules were down-regulated in sepsis compared to control (Fig. 2A). Modules M14, M63, and M64 were also significantly changed in non-survived patients compared with survived patients (Fig. 2B). The three modules could perfectly discriminate samples of sepsis and control in both datasets GSE185263 and GSE65682 (Fig. 3). Other modules including M2, M33, M35, M57, and M63 had an AUC larger than 0.9.
Modules were differentially expressed in SIRS, sepsis, septic shock, and ARDS
Severe sepsis or septic shock is characterized by an excessive inflammatory response to infectious pathogens [4]. We analyzed the dataset GSE63042 which contains transcriptome data for SIRS, sepsis, septic shock, and sepsis death. M33 was upregulated in severe sepsis, sepsis shock, and sepsis death compared to SIRS. M42, M52 was downregulated in severe sepsis and sepsis death compared to SIRS. M57 was upregulated in sepsis death compared to severe sepsis. M59 was upregulated in sepsis death compared to SIRS, uncomplicated sepsis, and septic shock. M63 was upregulated in severe sepsis, septic shock, and sepsis death compared to SIRS (Fig. 4).
ARDS is a devastating complication of severe sepsis, which results in a high mortality rate. We analyzed the dataset GSE32707 and found that M27, M49, M59, and M63 were upregulated and M30, M35, M39, M40, M42, M60, M66, and M70 were downregulated in day 0 ARDS compared to normal control.
Modules were associated with the survival of critically ill patients with sepsis
To check if the identified modules are associated with survival, we performed survival analysis in the dataset GSE65682 of critically ill patients with sepsis. Modules M17, M49, M52, and M59 were found to be associated with patient survival after adjusting for endotype class and age (Fig. 5). Hub genes of these modules were also associated with patient survival (Fig. 5). In sepsis dataset GSE185263, we found that M52 (R =-0.34, P =8E-11) and M60 (R =0.35, P =7E-12) had the highest absolute correlation value with Sequential Organ Failure Analysis (SOFA) score. M60 is enriched with non-coding genes. Survival analysis for all the modules can be explored at the sepsis gene co-expression database (https://liqs.shinyapps.io/sepsis/).
Module expression was changed in severe A (H1N1) pandemic influenza
Critically ill patients with sepsis caused directly by influenza viruses, or by influenza-induced secondary bacterial infections are increasing worldwide. Therefore, we performed differential analysis for GSE27131 to examine if these modules are associated with disease development in severe A (H1N1) pandemic influenza. ME59 were found to be strongly up-regulated in day 6 compared to day 0 (P =0.003), but not changed in day 0 compared to control. Thus, M59 was associated with severe influenza development. M59 hub gene CEACAM6 has been reported as a protein receptor for the influenza A virus [24]. Interestingly, CEACAM6-high airway neutrophils and epithelial cells are a feature of severe asthma [25]. In our analysis, the module was associated with neutrophil activation (P =5E-22), indicating CEACAM6 as a potential biomarker for severe influenza that requires mechanical ventilator support. The finding was confirmed in another dataset GSE32707, in which both day 0 and day 7 sepsis patients with acute respiratory distress syndrome (ARDS) had stronger M59 expression compared to control (P <0.0001) [9]. In M59, many genes were co-expressed with CEACAM6 (Fig. 6A), such as CEACAM8 which can also separate patients with long and short survival (Fig. 6B).
Modules can separate patients with severe community-acquired pneumonia
Diagnosis of severe pneumonia remains challenging because of a lack of correlation between the clinical status and etiology [26]. Causes of severe respiratory failure include bacterial pneumonia, influenza pneumonia, mixed bacterial and influenza A pneumonia, and non-infective SIRS. We used dataset GSE40012 to test if the modules can separate pneumonia patients with different etiologies. For each disease, we selected the top one or two modules to describe here. For example, M56 can separate bacterial pneumonia from other pneumonia (Fig. 7A). M17, M27, and M59 can separate influenza pneumonia from other pneumonia (Fig. 7BC). M69 can separate mixed-type pneumonia from other pneumonia (Fig. 7D). M67 can separate SIRS from other pneumonia (Fig. 7E). M14 can separate health control from pneumonia (Fig. 7F).
Candidate drugs associated with sepsis
Four modules M17, M49, M52, and M59 were associated with patient survival. We searched the genes of these modules in DSigDB and found that testosterone (P =3E-145) and phytoestrogens (P =2E-55) signatures were enriched with M17 genes. Ibuprofen (P =6E-4) and urea (P =8E-4) signatures were enriched with M49 genes. Dichlorvos (P =1E-6) and ZIRAM (P =1E-6) signatures were enriched with M52 genes. Potassium persulfate (P =2E-5) and adenylyl sulfate (P =4E-4) signatures were enriched with M59 genes. GWAS Catalog analysis revealed that M49 and M59 were associated with monocyte percentage of white cells (P =0.004) and vitamin B12 levels (P =0.007) respectively. We also performed docking analysis to hub gene of M49 AQP9 in complex with urea (Fig. 8A). The interaction may occur at sites Ala214 and Asn216 (Fig. 8B)..
A user-friendly database for sepsis gene co-expression was constructed
Although sepsis transcriptomes have been analyzed in several publications, there are still few web tools available for researchers. Thus, we for the first time constructed an easy-to-use web tool that provides informative data about the above sepsis gene co-expression network. The modules of the database includes module gene list, module hub gene list, module network visualization, gene function investigation, and differential module analysis (Fig. 9). The database is available at https://liqs.shinyapps.io/sepsis/. The database will provides a useful tool for researchers to generate testable hypothesis.