Screening of Important Factors in the Early Sepsis Stage Based on the Evaluation of ssGSEA Algorithm and ceRNA Regulatory Network.

BACKGROUND: Sepsis is a dysregulated host response to pathogens. Delay in sepsis diagnosis has become a primary cause of patient death. This study determines some factors to prevent septic shock in its early stage, contributing to the early treatment of sepsis.
METHODS: The sequencing data (RNA- and miRNA-sequencing) of patients with septic shock were obtained from the NCBI GEO database. After re-annotation, we obtained lncRNAs, miRNA, and mRNA information. Then, we evaluated the immune characteristics of the sample based on the ssGSEA algorithm. We used the WGCNA algorithm to obtain genes significantly related to immunity and screen for important related factors by constructing a ceRNA regulatory network. RESULT: After re-annotation, we obtained 1708 lncRNAs, 129 miRNAs, and 17 326 mRNAs. Also, through the ssGSEA algorithm, we obtained 5 important immune cells. Finally, we constructed a ceRNA regulation network associated with SS pathways.
CONCLUSION: We identified 5 immune cells with significant changes in the early stage of septic shock. We also constructed a ceRNA network, which will help us explore the pathogenesis of septic shock.

Introduction

Sepsis is a dysregulated host response to pathogens. It causes circulatory, cellular, or metabolic circulatory abnormalities and ultimately leads to life-threatening organ dysfunction. Septic shock (SS) is the most severe complication of sepsis and is also the primary cause of mortality in intensive care units (ICUs) worldwide. Early and appropriate regulations are crucial for patient recovery; therefore, it is necessary to explore SS-related factors to help early diagnosis, monitor, and prevent SS development. The pathogenesis of SS is complicated. Cells of immune systems contribute to its occurrence and development. After septic shock, T and B lymphocytes are altered.[3,4] Therefore, the proportion of circulating and regulatory cells will increase. Additionally, neutrophils play a critical role during septic shock. Early higher neutrophils counts relate to the increase in sepsis severity. Therefore, exploring changes in the type and proportion of immune cells in the early samples of SS will be helpful for rapid diagnosis.

Recently, long non-coding RNAs (lncRNAs) have become a hot spot in disease research. The mutual regulation mode between lncRNA, miRNA, and its downstream target genes is closely related to the occurrence and development of diseases. For example, Manetti et al found at least 77 miRNAs involved in cardiac inflammation and dysfunction during sepsis. Additionally, lncRNA can regulate miRNA (an important factor of post-transcriptional regulation activity) through sponge adsorption. This lncRNA is called competitive endogenous RNA (ceRNA). LncRNA (a ceRNA) can competitively bind to miRNA, thereby regulating the protein level of coding genes and the biological behavior of cells.[9,10] Therefore, an in-depth understanding of ceRNA-based regulatory mechanisms will help us better understand the pathogenesis of SS and screen important genes.

This research screened out RNAs whose expression levels changed significantly at 24 and 48 hours after SS and constructed a ceRNA regulatory network. Additionally, by gene set enrichment analysis (GSEA), we screened out immune cells with evident differences in different SS periods.

Material and Methods

Data collection

We downloaded the data set, GSE57065, from the NCBI GEO database. There were 82 samples related to SS in GSE57065, and these 82 blood samples were enrolled from 28 ICU patients at the onset of sepsis shock 0, 24, and 48 hours. The sample detection platform was GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array. According to the time point, we divided the sample into 2 comparison groups: 24 hours versus 0 hour and 48 hours versus 0 hour. Then, we downloaded the detailed annotation information of the GPL570 Affymetrix Human Genome U133 Plus 2.0 Array platform from the Ensembl genome browser 96 database (http://asia.ensembl.org/index.html), including the probe, gene Symbol, RNA type, and other information. Finally, we re-annotated the expression data in GSE57065 and obtained the lncRNA, miRNA, and mRNA expression levels based on these annotations.

Screening of significantly differentially expressed genes

After re-annotation, we used the limma package of the R software (https://bioconductor.org/packages/release/bioc/html/limma.html) to screen for significantly differentially expressed RNAs (DERs, including lncRNA, miRNA, and mRNA) in the 24 hours versus 0 hour and 48 hours versus 0 hour comparison groups. FDR < 0.05 and |log2FC| > 0.263 were considered statistically significant. Heatmaps were generated using pheatmap package (https://cran.r-project.org/web/packages/pheatmap/index.html), and volcano plots were also conducted in R software.[12,13] Finally, we retained the overlapping DERs in the 2 comparison groups. These DERs were genes whose expression levels continue to change significantly during the process of 24 to 48 h in the early stage of SS.

Immune-infiltrating cell analysis

Here, we downloaded the immunologic signature gene sets from GSEA database (http://software.broadinstitute.org/gsea/index.jsp). Then we use the ssGSEA algorithm in the GSVA package of the R software (http://www.bioconductor.org/packages/release/bioc/html/GSVA.html) to evaluate the type of immune infiltration of the samples in GSE57065. Finally, we compared the proportion of each immune cell in the 0, 24, and 48 hours samples and kept the immune cell types that were significantly different from the 0 hour time for subsequent analysis.

Screening disease progression and immunity-related modules by WGCNA algorithm

For all genes with expression detected in GSE57065, we used the WGCNA R package (https://cran.r-project.org/web/packages/WGCNA/index.html) to screen for significantly stable gene modules that are related to the sample stage and immune cells. Screening thresholds: contains at least 100 genes, cutHeight = 0.995.

Then, we mapped the overlapping genes screened in step 2.2 to each WGCNA module. Finally, we calculated the significant enrichment parameter-fold enrichment and enrichment significance P-value through the hypergeometric algorithm. f(k,N,M,n) = C(k,M)*C(nk,NM)/C(n,N).

Where, N represents all genes involved in the analysis of WGCNA; M represents the number of genes in each module obtained by WGCNA analysis; n represents the number of DERs obtained in step 2.2; and k represents the number of DERs in the intersection mapped to the corresponding module.

In this research, we selected P < .05 and fold enrichment > 1 as the thresholds. We took the significantly enriched genes in the target module as the object of further analysis and research.

Construction of the ceRNA network

We used DIANA-LncBase v.2 (http://www.microrna.gr/LncBase) to predict the lncRNA–miRNA interactions based on DELncRNAs obtained from step 2.2. StarBase v.2.0 (http://starbase.sysu.edu.cn/) was used to search for target genes for the DEmiRNA. Then, we mapped the genes that were significantly enriched in WGCNA modules to the target genes. Additionally, we used these genes to construct interactions. Ultimately, we integrated these interactions and constructed a lncRNA-related ceRNA network with Cytoscape v.3.6.1 based on how lncRNAs can affect the function of miRNAs and act as miRNA sponges to regulate mRNA expression.

We performed a functional analysis of mRNA in the ceRNA network using DAVID v.6.8 (annotation, visualization, and integrated discovery database, http://david.abcc.ncifcrf.gov/).[20,21]

Screening disease-related pathway

We used “septic shock” as the keyword to search for the SS-related Kyoto Encyclopedia for Genes and Genomes (KEGG) signal pathways in the Comparative Toxicogenomics Database (CTD) (http://ctd.mdibl.org/). Compared with the KEGG signal pathways obtained from step 2.5, we kept the overlapped pathways and constructed a regulated network associated with SS.

All data-related scripts were provided in Supplemental Material 1.

Result

After re-annotation, we got 1708 lncRNAs, 129 miRNAs, and 17 326 mRNAs. Then, 874 and 1546 DERs were selected by Limma in the 24 hours versus 0 hour and 48 hours versus 0 hour comparison groups, respectively. The volcano map (Figure 1 left) and heatmap (Figure 1 right) of the lncRNAs, miRNAs, and mRNAs showed that these DERs could separate different samples based on time. By comparing the DERs screened in the 2 comparison groups, 644 overlapped DERs were retained, including 26 lncRNAs, 7 miRNAs, and 611 mRNAs. These DERs were RNAs whose expression levels were significantly different in the 24 and 48 hours samples.

Figure 1.

DEGs identified in 24 h versus 0 h (A) and 48 h versus 0 h (B). Left: Volcano map. Red dots: upregulated genes; Blue dots: downregulated genes. The horizontal red dotted line indicates false discovery rate (FDR) = 0.05, and the 2 vertical red dotted lines indicate |log2 FC| (fold change) = 0.263. Right: Heatmap.

Based on the gene expression data, using the ssGSEA algorithm, we obtained 28 immune cell ratios. First, we compared the differences in the proportion of various immune cells at 0, 24, and 48 h. Then, we kept the immune cell types that were significantly different at 24 h or 48 h time points compared with the time at 0 h, and finally got 5 immune cells with significant differences (Figure 2). The 5 immune cells were activated CD4 T cell, T follicular helper cell, Regulatory T cell, activated B cell, and neutrophils.

Figure 2.

Distribution diagram of immune cell types with significant differences between the 0, 24, and 48 h time points.

*P < .05. **.005 < P < .05. ***P < 0.005.

Module detection

First, we assume that the gene network is subject to scale-free distribution. Then, we selected a soft threshold (power) to make the constructed network a scale-free network. In this research, 9 was set as the soft threshold to meet the selected criteria of power value (Figure 3A). Based on the standard of a dynamic cut tree, 100 genes were set as the least gene number of each gene networks and 0.995 as the cut height, respectively. Finally, 13 modules were screened out.

Figure 3.

(A) Determination of soft-threshold power in the WGCNA. (Left) Analysis of the scale-free index for various soft-threshold powers. (Right) Analysis of the mean connectivity for various soft-threshold powers. (B) Dendrogram of different modules. Each color represents a different module.

Based on results of step 3.2, the correlation between modules, immune-infiltrating cells, and time was calculated, and modules related to different immune-infiltrating cells and time were screened (Figure 4).

Figure 4.

Identification of modules associated with immune cells and time.

Based on the hypergeometric enrichment algorithm described in the method, we mapped 611 overlapped DERs to WGCNA modules. Therefore, 277 DERs are located in WGCNA modules (As shown in Table 1). Furthermore, they were significantly enriched in black and magenta modules, containing 73 and 43 genes, respectively.

Table 1.

Modules.

ID	Color	Module size	# Overlapped DEGs	Enrichment information
				Enrichment fold [95%CI]	Phyper
Module 1	Black	205	73	6.943 [5.104-9.377]	2.20E-16
Module 2	Blue	414	12	0.566 [0.286-1.015]	5.75E-02
Module 3	Brown	341	14	0.801 [0.428-1.386]	5.23E-01
Module 4	Green	223	6	0.525 [0.189-1.176]	1.53E-01
Module 5	Greenyellow	118	1	0.165 [0.00414-0.947]	4.61E-02
Module 6	Gray	2388	90	0.735 [0.570-0.941]	1.25E-02
Module 7	Magenta	169	43	4.496 [3.390-7.140]	1.32E-14
Module 8	Pink	176	9	0.998 [0.444-1.966]	9.99E-01
Module 9	Purple	161	2	0.242 [0.0289-0.898]	2.43E-02
Module 10	Red	217	1	0.0899 [0.00226-0.510]	4.80E-04
Module 11	Tan	113	1	0.173 [0.00432-0.990]	4.41E-02
Module 12	Turquoise	560	25	0.871 [0.549-1.327]	6.12E-01
Module 13	Yellow	320	–	–	–

The modules in bold were the modules with significantly enriched DERs.

Prediction of lncRNA-miRNA-mRNA interactions

First, we used the experimental module DIANA-LncBase v.2 (http://www.microrna.gr/LncBase) to predict the lncRNA–miRNA interactions based on the overlapped DELncRNAs. Then, T starBase v.2.0 (http://starbase.sysu.edu.cn/) was used to predict the interactions between DEmiRNAs and DEmRNAs. Next, we mapped DERNAs enriched in black and magenta modules to regulate target genes and predict miRNA–mRNA interactions. Finally, we constructed a lncRNA-related ceRNA network using Cytoscape v.3.6.1 (Figure 5). It shows which lncRNAs can affect the function of miRNAs and how they regulate mRNA expression.

Figure 5.

ceRNA network. Squares, triangles, and circles represent DElncRNA, DEmiRNA, and DEmRNA, which are significantly enriched in modules.

Gene function analysis

We used DAVID v.6.8 for KEGG signaling pathway enrichment annotation analysis on mRNAs in the ceRNA regulatory network, and 1 of 9 KEGG signal pathways was obtained (Table 2).

Table 2.

KEGG signaling pathway with significant mRNA correlation in the ceRNA regulatory network.

Term	Count	P-Value
hsa00410: beta-Alanine metabolism	4	4.13E-05
hsa00360: phenylalanine metabolism	2	7.86E-03
hsa00340: histidine metabolism	2	1.01E-02
*hsa01100: metabolic pathways	10	1.16E-02
hsa00350: tyrosine metabolism	2	1.55E-02
hsa00260: glycine, serine, and threonine metabolism	2	1.71E-02
hsa00071: fatty acid degradation	2	1.83E-02
hsa00280: valine, leucine, and isoleucine degradation	2	2.03E-02
*hsa04060: cytokine-cytokine receptor interaction	3	3.26E-02

Screening disease-related pathway

We used “septic shock” as a keyword in the CTD database and searched 75 KEGG signal pathways directly related to SS. Compared with the pathways obtained before, Metabolic and cytokine–cytokine receptor interaction signaling pathways were screened. Then, we separately extracted the parts of the ceRNA regulatory network directly related to these 2 disease pathways (Figure 6). In these 2 pathways, has-miR-4668, has-miR-6847, has-miR-601, and has-miR-6281 contribute to SS pathogenesis. GABPB1-AS1 can regulate most of the miRNAs in the ceRNA network.

Figure 6.

ceRNA network directly associates with SS. Squares, triangles, and circles represent DElncRNA, DEmiRNA, and DEmRNA, which significantly enriched in metabolic pathways and cytokine–cytokine receptor interaction signaling pathways, respectively.

Discussion

Sepsis is mostly caused by the dysregulated host response to infection. SS is a life-threatening condition that could complicate sepsis. Septic shock occurs in more than 230 000 US patients yearly. Septic shock requires prompt diagnosis, especially for severe infection. In this study, we downloaded SS expression profile data from the NCBI GEO database and re-annotated them. We classified the samples according to time (0, 24, and 48 hours) and filtered out 644 DERs to find new clues for SS diagnosis.

Septic shock results from the dysregulation of the innate immune response. The immune microenvironment significantly affects the diagnosis, survival outcome, and clinical treatment sensitivity of the disease. Additionally, the immune microenvironment is composed of immune cells, immune-related pathways, and cytokines secreted by immune cells; therefore, different immune cell-infiltration patterns may provide a favorable or unfavorable environment for disease development. Our research screened 5 important immune cells: activated CD4 T cells, T follicular helper cells, regulatory T cells, activated B cells, and neutrophils. The ratio of T follicular helper cells and regulatory T cells gradually decreased over time. The ratio of activated B cells and neutrophils gradually increased over time. Interestingly, the proportion of activated CD4 T cells will first decrease and then increase. At the beginning of sepsis, the number of neutrophils will increase significantly.[24,25] The adaptive immune system is composed of T and B lymphocytes ; therefore, we can diagnose SS by detecting changes in the proportion of these 5 immune cells.

To identify the biological key genes related to SS, we applied WGCNA to select the DEGs. We found that the DEGs were significantly enriched in black and magenta modules. Based on the DEGs in these modules, we constructed a ceRNA network. Next, KEGG pathway analyses were performed. Compared with pathways selected from CTD databases, metabolic pathways and cytokine–cytokine receptor interaction signaling pathways are present. These 2 pathways contain many important genes, such as ABAT, IL1RAP, OSM, and PIGX. IL1RAP plays an important role in blood disease and is emerging as a novel therapeutic target.[26,27] ABAT and PIGX could suppress the growth of some tumor cells.[28-30] The expression of OSM will be upregulated in patients with sepsis. In ceRNA regulation network based on these DERs, has-miR-4668 plays an important central regulatory role. Ten mRNAs in the network were regulated by hsa-miR-4668; hsa-miR-4668 was also associated with ten lncRNAs. A previous study showed that MIR4668 is induced by B. pseudomallei infection, and inhibits autophagy, conferring bacterial subversion of host innate immune pathways.

Conclusions

Septic shock is a clinical emergency and needs rapid diagnosis; however, the diagnosis of shock is complex and more difficult for severe infection. This review found some bio-markers that help rapid diagnosis of SS. The ratio of activated CD4 T cells, T follicular helper cells, regulatory T cells, activated B cells, and neutrophils will change significantly after septic shock. This provides us with new ideas for the rapid diagnosis of SS.

Click here for additional data file.

Supplemental material, sj-txt-1-evb-10.1177_11769343211058463 for Screening of Important Factors in the Early Sepsis Stage Based on the Evaluation of ssGSEA Algorithm and ceRNA Regulatory Network by Liou Huang, Chunrong Wu, Dan Xu, Yuhui Cui and Jianguo Tang in Evolutionary Bioinformatics