Identification of hub genes and prediction of the ceRNA network in adult sepsis

Abstract

Background: Sepsis refers to a dysregulated host immune response to infection. It carries a high risk of morbidity and mortality, and its pathogenesis has yet to be fully elucidated. The main aim of this study was to identify prognostic hub genes for sepsis and to predict a competitive endogenous RNA (ceRNA) network that regulates the hub genes.

Methods: Six transcriptome datasets from the peripheral blood of septic patients were retrieved from the Gene Expression Omnibus (GEO) database. The robust rank aggregation (RRA) method was used to screen differentially expressed genes (DEGs) across these datasets. A comprehensive bioinformatics investigation was conducted, encompassing Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses using the "clusterProfiler" package in R, as well as gene set enrichment analysis (GSEA) to further elucidate the biological functions and pathways associated with the DEGs. Weighted gene co-expression network analysis (WGCNA) was performed to identify a module significantly associated with sepsis. Integration of this module with protein-protein interaction (PPI) network analysis facilitated the identification of five hub genes. These hub genes were subsequently validated using an independent dataset and reverse transcription-quantitative polymerase chain reaction (RT-qPCR) analysis of peripheral blood samples from septic patients. The prognostic values of these hub genes were assessed via receiver operating characteristic (ROC) curve analysis. Finally, a ceRNA network regulating the prognostic hub genes was constructed by integrating data from a literature review as well as five online databases.

Results: RRA analysis identified 164 DEGs across six training cohorts. Bioinformatics analyses revealed concurrent hyperinflammation and immunosuppression in sepsis patients. Five hub genes were identified via WGCNA and PPI network analysis, and their differential expression was verified by the validation dataset (GSE28750) and RT-qPCR analysis in the peripheral blood of septic patients. ROC analysis confirmed four hub genes with prognostic value, and a ceRNA network was predicted to elucidate their regulatory mechanisms.

Conclusion: This study identified four hub genes (CLEC4D, GPR84, S100A12, and HK3) with significant prognostic value in sepsis and predicted a ceRNA network (NEAT1-hsa-miR-495-3p-ELF1) regulating their expression. The integrated analysis reconfirmed the concurrent presence of hyperinflammation and immunosuppression in hospitalized sepsis patients. These findings enhance the understanding of sepsis pathogenesis and identify potential therapeutic targets.

Keywords: Inflammation; Peripheral blood; Prognostic biomarker; RRA; WGCNA.

Figure 1. Flowchart of the study.

Figure 2. Identification of DEGs by RRA analysis.

(A) Volcano plots showing DEGs between septic patients and healthy controls across six Gene Expression Omnibus datasets. Red points indicate upregulated genes in sepsis, while blue points indicate down-regulated genes in sepsis. Gray points indicate genes with no significant difference. (B) Heatmap showing the top 10 up- and down-regulated DEGs in RRA. Red indicates highly expressed genes in sepsis, while green indicates poorly expressed genes in sepsis. (C) The enriched biological functions in Gene Ontology analysis demonstrated the underlying significance of up- and down-regulated DEGs across biological process, cellular component, and molecular function categories. (D) Top 10 pathways enriched in up- and down-regulated DEGs, as according to Kyoto Encyclopedia of Genes and Genomes. DEGs, differentially expressed genes; RRA, robust rank aggregation.

Figure 3. Gene set enrichment analysis annotation of statistically changed genes in robust rank aggregation analysis.

(A) Ridgeline plot displaying the distribution of gene expression in the top six most significant gene sets. (B) Gene sets enriched in neutrophil degranulation (NES = 3.104, P = 1 × 10⁻¹⁰, q = 2.005 × 10⁻⁸), innate immune system (NES = 2.604, P = 1  ×  10⁻¹⁰, q = 2.005 × 10⁻⁸), second messenger molecules (NES = −2.969, P = 3.879 ×  10⁻⁸, q = 3.185 × 10⁻⁶), co-stimulation by CD28 family (NES = −2.734, P = 2.002 ×  10⁻⁶, q = 2.027 ×  10⁻⁴), diseases of the immune system (NES = 2.300, P = 5.65 × 10⁻⁶, q = 4.253 × 10⁻⁴), metabolism of RNA (NES = 1.897, P = 1.479 × 10⁻⁵, q = 0.001). NES, normalized enrichment score.

Figure 4. Estimation of the infiltration of immune cells in the six training datasets using the CIBERSORT analysis.

(A–F) Stacked bar charts revealed the percentage of immune cells in six datasets (A) GSE137340, (B) GSE69063, (C) GSE69528, (D) GSE54514, (E) GSE57065 and (F) GSE95233, and boxplots revealed the differential infiltration of immune cells between patients with sepsis and healthy controls in six datasets (adjusted P-values). ^∗P < 0.05, ^∗∗P < 0.01, ^∗∗∗P < 0.001.

Figure 5. Weighted gene co-expression network analysis in the GSE69528 dataset.

(A) Establishment of soft-thresholding power and its mean connectivity. (B) Hierarchical clustering tree showing 11 cluster modules. (C) Heatmap showing the correlation of each module with sepsis. Rows represent modules while columns represent clinical phenotypes. Red indicates a positive correlation, while blue indicates a negative correlation. (D) Scatterplot showing the association between gene significance and module membership in the blue module. (E) Key genes in the blue module were subjected to gene ontology enrichment analysis. (F) Key genes in the blue module were subjected to kyoto encyclopedia of genes and genomes enrichment analysis. (G) Venn diagram showing the obtained gene set.

Figure 6. PPI network construction and hub gene identification.

(A) A PPI network was constructed and visualized using Cytoscape. The nodes represent proteins, while the edges represent their interactions. (B) The plug-in Minimal Common Oncology Data Elements screened two subnetworks: Cluster 1 and Cluster 2. (C and D) The plug-in cytoHubba- Maximal Clique Centrality was performed to identify the top 10 genes in the network. (E) The plug-in iRegulon highlights the predicted transcription factor in blue, while target genes are shown in orange. (F) Analysis of differential gene expression in the GSE28750 validation dataset. ^∗∗∗P < 0.001. (G) The expression of hub genes was validated in peripheral blood samples from patients with sepsis and healthy controls using reverse transcription-quantitative PCR. ^∗∗P < 0.01, ^∗∗∗P < 0.001. (H) Receiver operating characteristic curve showing prognosis. PPI, protein–protein interaction.

Figure 7. Correlation between the prognostic hub genes and immune cell types as well as immune genes.

(A) Correlation analysis of the four prognostic hub genes. (B) Analysis of the correlation between the four prognostic hub genes and immune cells. (C) Analysis of the correlation between the four prognostic hub genes and immune genes. ^∗P < 0.05, ^∗∗P < 0.01.

Figure 8. Construction of a miRNA-mRNA network and a ceRNA network.

(A) ELF1 (blue hexagon) and 98 regulatory miRNAs (pink circle) formed the miRNA-mRNA network. (B) Construction of a ceRNA network where the green elliptical, pink round node and blue hexagonal nodes represent long noncoding RNA, miRNA and mRNA, respectively . miRNA, microRNA; ceRNA, competitive endogenous RNA.