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

Abstract

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.

Keywords

Sepsis pathogenicity genes regulator

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 |log₂FC| > 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

Screening of significantly differentially expressed genes

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 |log₂ FC| (fold change) = 0.263. Right: Heatmap.

Immune-infiltrating cell analysis

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.

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
ID	Color	Module size	# Overlapped DEGs	Enrichment fold [95%CI]	P_hyper
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.

Construction of the ceRNA network

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.

Supplemental Material

sj-txt-1-evb-10.1177_11769343211058463 – Supplemental material for Screening of Important Factors in the Early Sepsis Stage Based on the Evaluation of ssGSEA Algorithm and ceRNA Regulatory Network

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

Footnotes

Funding:

The author(s) received no financial support for the research, authorship, and/or publication of this article.

Declaration of Conflicting Interests:

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Authors’ Contributions

Jianguo Tang and Liou Huang were responsible for the conception and design of the research, and drafting the manuscript. Chunrong Wu and Dan Xu performed the data acquisition. Liou Huang and Yuhui Cui performed the data analysis and interpretation. All authors have read and approved the manuscript.

ORCID iD

Jianguo Tang

Supplemental material

Supplemental material for this article is available online.

References

Esposito

De Simone

Boccia

De Caro

Pagliano

. Sepsis and septic shock: New definitions, new diagnostic and therapeutic approaches. J Glob Antimicrob Resist. 2017;10:204-212.

Annane

Bellissant

Cavaillon

. Septic shock. Lancet. 2005;365:63-78.

Gustave

Gossez

Demaret

, et al. Septic shock shapes B cell response toward an exhausted-like/immunoregulatory profile in patients. J Immunol. 2018;200:2418-2425.

Rimmelé

Payen

Cantaluppi

, et al. Immune cell phenotype and function in sepsis. Shock. 2016;45:282-291.

Stiel

Meziani

Helms

. Neutrophil activation during septic shock. Shock. 2018;49:371-384.

Park

Kim

Choe

, et al. Neutrophils disturb pulmonary microcirculation in sepsis-induced acute lung injury. Eur Respir J. 2019;53:1800786.

Manetti

Maiese

Paolo

, et al. MicroRNAs and sepsis-induced cardiac dysfunction: a systematic review. Int J Mol Sci. 2020;22:321.

Wang

Liu

, et al. CREB up-regulates long non-coding RNA, HULC expression through interaction with microRNA-372 in liver cancer. Nucleic Acids Res. 2010;38:5366-5383.

Liu

Sun

Nie

, et al. Lnc RNA HOTAIR functions as a competing endogenous RNA to regulate HER2 expression by sponging miR-331-3p in gastric cancer. Mol Cancer. 2014;13:92.

10.

Zhou

Gao

Wang

Zhang

Liu

Duan

. Linc-RNA-RoR acts as a “sponge” against mediation of the differentiation of endometrial cancer stem cells by microRNA-145. Gynecol Oncol. 2014;133:333-339.

11.

Ritchie

Phipson

, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43:e47.

12.

Wang

Cao

, et al. RNA-seq analyses of multiple meristems of soybean: novel and alternative transcripts, evolutionary and functional implications. BMC Plant Biol. 2014;14:169.

13.

Eisen

Spellman

Brown

Botstein

. Cluster analysis and display of genome-wide expression patterns. Proc Natl Acad Sci USA. 1998;95:14863-14868.

14.

Zhang

Kang

, et al. Tumor-infiltrating immune cells act as a marker for prognosis in colorectal cancer. Front Immunol. 2019;10:2368.

15.

Langfelder

Horvath

. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.

16.

Cao

Zhang

. A Bayesian extension of the hypergeometric test for functional enrichment analysis. Biometrics. 2014;70:84-94.

17.

Paraskevopoulou

Vlachos

Karagkouni

, et al. DIANA-LncBase v2: indexing microRNA targets on non-coding transcripts. Nucleic Acids Res. 2016;44:D231-D238.

18.

Liu

Tan

Liang

, et al. A reduction in reactive oxygen species contributes to dihydromyricetin-induced apoptosis in human hepatocellular carcinoma cells. Sci Rep. 2014;4:7041.

19.

Shannon

Markiel

Ozier

, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498-2504.

20.

Huang

Sherman

Lempicki

. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4:44-57.

21.

Huang

Sherman

Lempicki

. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009;37:1-13.

22.

Davis

Grondin

Johnson

, et al. Comparative Toxicogenomics Database (CTD): update 2021. Nucleic Acids Res. 2021;49:D1138-d1143.

23.

Huet

Chin-Dusting

. Septic shock: desperately seeking treatment. Clin Sci. 2014;126:31-39.

24.

Heffernan

Monaghan

Thakkar

Machan

Cioffi

Ayala

. Failure to normalize lymphopenia following trauma is associated with increased mortality, independent of the leukocytosis pattern. Crit Care. 2012;16:R12.

25.

Hoser

Skirecki

Złotorowicz

Zielińska-Borkowska

Kawiak

. Absolute counts of peripheral blood leukocyte subpopulations in intraabdominal sepsis and pneumonia-derived sepsis: a pilot study. Folia Histochem Cytobiol. 2012;50:420-426.

26.

Mitchell

Barreyro

Todorova

, et al. IL1RAP potentiates multiple oncogenic signaling pathways in AML. J Exp Med. 2018;215:1709-1727.

27.

Warda

Larosa

Neto

Rocha

, et al. CML hematopoietic stem cells expressing IL1RAP can be targeted by chimeric antigen receptor-engineered T cells. Cancer Res. 2019;79:663-675.

28.

Shen

Wang

Ying

McAlinden

O’Keefe

. Inhibition of 4-aminobutyrate aminotransferase protects against injury-induced osteoarthritis in mice. JCI Insight. 2019;4:e128568.

29.

Chen

Cao

Liao

, et al. Loss of ABAT-mediated GABAergic system promotes basal-like breast cancer progression by activating Ca(2+)-NFAT1 axis. Theranostics. 2019;9:34-47.

30.

Nakakido

Tamura

Chung

, et al. Phosphatidylinositol glycan anchor biosynthesis, class X containing complex promotes cancer cell proliferation through suppression of EHD2 and ZIC1, putative tumor suppressors. Int J Oncol. 2016;49:868-876.

31.

Hua

Zhang

Gong

, et al. Dihydromyricetin improves endothelial dysfunction in diabetic mice via oxidative stress inhibition in a SIRT3-dependent manner. Int J Mol Sci. 2020;21:6699.

32.

Fang

Zhu

, et al. Burkholderia pseudomallei survival in lung epithelial cells benefits from miRNA-mediated suppression of ATG10. Autophagy. 2015;11:1293-1307.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB