Comprehensive Analysis of Molecular Subtypes and Hub Genes of Sepsis by Gene Expression Profiles

Abstract

Background: Sepsis is a systemic inflammatory response syndrome (SIRS) with heterogeneity of clinical symptoms. Studies further exploring the molecular subtypes of sepsis and elucidating its probable mechanisms are urgently needed.

Methods: Microarray datasets of peripheral blood in sepsis were downloaded from the Gene Expression Omnibus (GEO) database, and differentially expressed genes (DEGs) were identified. Weighted gene co-expression network analysis (WGCNA) analysis was conducted to screen key module genes. Consensus clustering analysis was carried out to identify distinct sepsis molecular subtypes. Subtype-specific pathways were explored using gene set variation analysis (GSVA). Afterward, we intersected subtype-related, dramatically expressed and module-specific genes to screen consensus DEGs (co-DEGs). Enrichment analysis was carried out to identify key pathways. The least absolute shrinkage and selection operator (LASSO) regression analysis was used for screen potential diagnostic biomarkers.

Results: Patients with sepsis were classified into three clusters. GSVA showed these DEGs among different clusters in sepsis were assigned to metabolism, oxidative phosphorylation, autophagy regulation, and VEGF pathways, etc. In addition, we identified 40 co-DEGs and several dysregulated pathways. A diagnostic model with 25-gene signature was proven to be of high value for the diagnosis of sepsis. Genes in the diagnostic model with AUC values more than 0.95 in external datasets were screened as key genes for the diagnosis of sepsis. Finally, ANKRD22, GPR84, GYG1, BLOC1S1, CARD11, NOG, and LRG1 were recognized as critical genes associated with sepsis molecular subtypes.

Conclusion: There are remarkable differences in and enriched pathways among different molecular subgroups of sepsis, which may be the key factors leading to heterogeneity of clinical symptoms and prognosis in patients with sepsis. Our current study provides novel diagnostic and therapeutic biomarkers for sepsis molecular subtypes.

Introduction

Sepsis, a systemic inflammatory response syndrome (SIRS) caused by various infectious processes, is one of the common diseases that leads to the death of hospitalized patients in the ICU (Fleischmann et al., 2016; Novosad et al., 2016). At present, we lack effective strategies for early diagnosis and treatment of sepsis due to the heterogeneity of pathogenesis and clinical symptoms in patients with sepsis (Davenport et al., 2016; Seymour et al., 2019). In addition, with the progression of sepsis and the persistence of systemic inflammation, patients with severe sepsis are commonly accompanied by multiple organ dysfunction syndromes (MODS), hypoperfusion, or hypotension, which brings a great challenge to sepsis treatment (Rello et al., 2017). Therefore, a thorough understanding of the distinct molecular subtypes of sepsis is needed to improve the prognosis of sepsis patients with different subtypes.

The release of a large number of cytokines and inflammatory mediators can result in dysregulated immune responses, which may be the decisive factor affecting the prognosis of sepsis patients (Markwart et al., 2015). In addition, the activation of abnormal genes in sepsis patients may play a critical role in the progression of the disease. In recent years, a variety of biomarkers, including serum cytokine/chemokine, cellular receptor, coagulation factors, vascular endothelial damage factors, and acute inflammatory factors have been implemented in the diagnosis and prognosis of sepsis (Pierrakos and Vincent, 2010; Sandquist and Wong, 2014). Nevertheless, because of the complexity of the pathogenesis of sepsis, the specificity and sensitivity of these biomarkers in disease diagnosis and prognosis are significantly lower than expected. Further exploring more potent biomarkers for the early diagnosis and treatment of sepsis has become urgent.

Bioinformatics analysis at the molecular biology level has been extensively applied for clinical practice by screening and predicting potential key pathways and biomarkers (Banwait and Bastola, 2015; Cheng et al., 2021). Few studies to date have explored the significance of identified molecular subtypes in the early diagnosis and treatment of diseases including cancer, respiratory diseases, and myocardial infarction (Hu et al., 2021; Li et al., 2021; Rao et al., 2021; Shi et al., 2021). Zhang et al. classified sepsis patients into three clusters based on m6A methylation regulatory genes (Zhang et al., 2020). However, whether the specific molecular subtypes could be determined based on whole genome sequencing data of sepsis patients is not yet thoroughly understood.

In this study, we classified sepsis into three molecular subtypes using unsupervised consensus clustering based on whole gene expression. Moreover, we identified consensus differentially expressed genes (co-DEGs) by intersecting the DEGs among three subtypes with differential genes screened by DEGs and WGCNA methods. Based on these results, we performed various analyses, including Gene set variation analysis (GSVA), Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), Reactome enrichment analysis, Pearson correlation analysis, and protein—protein interaction (PPI) analysis. In addition, we constructed a 25-gene-based diagnosis model using least absolute shrinkage and selection operator (LASSO) regression analysis and validated their expression levels and diagnostic values for sepsis. Eventually, we identified 7 distinct hub genes among the three molecular subtypes.

Materials

Data Acquisition

All GEO datasets were downloaded from Gene Expression Omnibus (GEO, www.ncbi.nlm.nih.gov/geo/) (Barrett et al., 2013). We selected 6 datasets (GSE154918, GSE54514, GSE9960, GSE69063, GSE25504, and GSE13904) related to sepsis for analysis. The whole gene expression profiles of peripheral blood were extracted for further analysis. The GEO datasets collected are exhibited in Table 1. The flowchart of the study was elucidated in Figure 1.

FIGURE 1

The GSE154918 dataset (GPL20301 platform) is composed of whole gene expression profiles of peripheral blood from 40 control and 24 sepsis samples, The GSE54514 dataset (GPL6947 platform) was composed of whole gene expression profiles of peripheral blood from 36 non-sepsis healthy control subjects and 127 sepsis patients, The GSE9960 dataset (GPL570 platform) was composed of whole gene expression profiles of peripheral blood from 16 control and 54 sepsis samples, The GSE69063 dataset (GPL20301 platform) was composed of whole gene expression profiles of peripheral blood from 33 control and 57 sepsis samples, The GSE25504 dataset (GPL570 platform) was composed of whole gene expression profiles of peripheral blood from 37 control and 26 sepsis samples, The GSE13904 dataset (GPL570 platform) was composed of whole gene expression profiles of peripheral blood from 18 control and 52 sepsis samples.

DEGs Analysis

R software’s “limma” package was applied for identifying the DEGs between the sepsis and control samples in the GSE154918 and GSE25504 datasets, respectively (Ritchie et al., 2015). DEGs with |log2 fold change (FC)| > 0.5 and p < 0.05 were defined as statistically significant. Volcano plots and heatmaps of the identified DEGs were visualized using the “ggplot2” and “heatmap” R packages.

Construction of the Co-Expression Network by WGCNA

R software’s “WGCNA” package was used to construct the co-expression network of the GSE154918 and GSE25504 datasets, respectively (Langfelder and Horvath, 2008). In brief, we explored the association among different pairs of genes and weighted them based on the expression levels of related genes in control and sepsis samples. Afterward, we transformed the adjacency matrix into a topological overlap matrix (TOM) to further verify the gene connectivity in the network. Next, the hierarchical clustering method was conducted to determine remarkably interconnected clusters (modules) according to gene connectivity and covariance coefficients. We selected the best power value and established the correlated modules. Finally, we performed Pearson correlation analysis and defined the highest three gene modules as statistically significant based on the correlation coefficient.

Unsupervised Consensus Clustering Analysis

The GSE9960, GSE13904, and GSE54514 datasets were log2 transformed. The robust multi-chip average (RMA) method was utilized to normalize gene expression, and the batch effect between different platforms was eliminated using the ComBat method based on the “sva” R package. Principal component analysis (PCA) was applied for evaluating the aggregation between different samples. Afterward, R software’s “ConsensusClusterPlus” package was carried out to perform unsupervised consensus analysis (Wilkerson and Hayes, 2010) among the 233 standardized sepsis patients. The number of re-samplings settled at 1,000, with each re-sampling containing 80% of the samples. The maximum number of clusters was set to 10, and the optimal k value was determined by the cumulative distribution function (CDF) index and the consensus matrix. Finally, t-SNE was carried out to validate the subtype assignments based on the gene expression profiles of the above sepsis patients.

The dataset was normalized and summarized using robust multi-chip average (RMA) implemented in the R package affy, and batch effects were corrected.

GSVA Analysis in Different Molecular Subtypes of Sepsis

R software’s “GSVA” and “GSEABase” packages were utilized to validate the performance of gene sets (c2. cp.kegg.v7.4. symbols) among sepsis patients with different molecular subtypes, thus identifying the enriched pathways in each subtype. Gene sets with adjusted p-value < 0.05 were defined as significantly enriched gene sets.

Identification and Correlation Analysis of Co-DEGs

The DEGs were screened using DEGs analysis and the WGCNA method in GSE154918 and GSE25504 datasets, and the specific genes were also identified in sepsis patients with distinct molecular subtypes. Eventually, a total of 48 core genes were identified by intersecting all the results. The “ggVennDiagram” package was used in generating the Venn diagrams of co-DEGs.

Spearman correlation analysis was performed to determine the correlation between core genes based on the gene expression profiles. The heatmap of the correlation coefficient among these hub genes was visualized using the “corrplot” R package, and the gene relationship network diagram with a correlation coefficient >0.9 was constructed using the “igraph” R package.

Enrichment Analysis

The Database for Annotation, Visualization, and Integrated Discovery (DAVID, https://david.ncifcrf.gov/summary.jsp), an online tool, was utilized to perform the GO enrichment analysis of co-DEGs (Huang da et al., 2009). Statistically significant GO terms (BP) with FDR <0.05 were screened and the results were visualized using the “GOplot” R package.

KEGG and Reactome enrichment analysis of these co-DEGs were conducted using the “clusterProfile” and “org.Hs.eg.db” R packages (Yu et al., 2012). The statistically significant enrichment pathways with adjusted p-value < 0.05 were defined and visualized using a bubble plot.

PPI Network Analysis

The STRING database (https://strin g-db.org/) was applied for constructing the PPI network of co-DEGs, and the core genes with a combined score of more than 0.4 were screened (Szklarczyk et al., 2015). The protein–protein interaction information was visualized utilizing the Cytoscape software (version 3.8.2).

LASSO Regression Analysis

LASSO is a regularization method with strong predictability, that is better than regression analysis when examining high-dimensional data (Bader and Hogue, 2003). A total of 303 samples including 70 health and 233 sepsis patients from GSE9960, GSE13904, and GSE54514 datasets were applied for LASSO regression analysis. As a training set, 70% of samples (N = 217, 58 controls, and 169 sepsis samples) were randomly selected. As a validation set, 30% of samples (N = 86, 22 controls, and 54 sepsis samples) were selected. The LASSO regression analysis was performed based on the expression profiles of co-DEGs using the “glmnet” R package. The predictive model constructed in the training set was verified in the testing set. In addition, receiver operating characteristic (ROC) curves were drawn in the validation set and train set using the “pROC” R package to assess the performance of the constructive model. AUC more than 0.75 was defined as a model with a high diagnostic value.

Validation of Hub Genes

GSE154918 and GSE69063 datasets were selected as the validation set to verify the diagnostic value of these genes identified by LASSO in sepsis. Genes with AUC >0.9 in both the GSE154918 and GSE69063 were defined as the ultimate hub genes. Subsequently, we verified the expression of these hub genes in GSE154918 and GSE69063 datasets and in distinct molecular subtypes of sepsis.

Results

Identification of DEGs and Construction of Co-Expression Network

The flowchart of our study is shown in Figure 1. To identify sepsis-related genes, we firstly analyzed the DEGs between sepsis and control samples in the GSE154918 dataset. A total of 1,671 up-regulated and 1,623 down-regulated genes were determined using the DEGs method (Figures 2A,B, Supplementary Figure S1A). Subsequently, we applied the WGCNA method to study the co-expression network in the GSE154918 dataset. When the best soft threshold power settled at 6, the mean connectivity was more effective (Figures 2C,D). A total of 13 co-expressed gene modules were created. among which the blue (4,036 genes, R = 0.87, p < 0.05), turquoise (4,618 genes, R = 0.83, p < 0.05), and yellow modules (1719 genes, R = 0.8, p < 0.05) had the most correlation with sepsis patients (Figure 2E, Supplementary Figure S1B). In addition, the blue (R = 0.85, p < 0.05), turquoise (R = 0.95, p < 0.05), and yellow (R = 0.92, p < 0.05) modules were significantly correlated with module-related genes (Supplementary Figures S1C–H). Meanwhile, we identified 945 up-regulated and 793 down-regulated genes using the DEGs method in the GSE25504 dataset (Figures 3A,B, Supplementary Figure S2A). Based on the WGCNA method, a relatively high mean connectively was maintained by setting the best soft threshold power to 8 (Figures 3C,D). We clustered a total of 14 co-expressed gene modules, among which the blue (2,174 genes, R = 0.79, p < 0.05), brown (1717 genes, R = -0.83, p < 0.05), and turquoise (2,478 genes, R = 0.8, p < 0.05) modules were closely correlated with sepsis patients (Figure 3E, Supplementary Figure S2B). Moreover, blue (R = 0.88 p < 0.05), brown (R = 0.93, p < 0.05), and turquoise (R = 0.9, p < 0.05) modules were also remarkably correlated with module-related genes (Supplementary Figures S2C–H).

FIGURE 2

FIGURE 3

Consensus Clustering Analysis for Sepsis

The GSE9960, GSE13904, and GSE54514 datasets were mixed into a combined dataset containing a total of 233 sepsis samples. These three datasets exhibited obvious separation before batch correction (Figure 4A). While the batch effect among these datasets from different platforms had been successfully eradicated after batch correction (Figure 4B). Then, we performed molecular subtypes analysis in a combined dataset based on the expression of genes via the “ConsensusClusterPlus” R package. Clustering results suggested that the classification was most reliable and stable when k = 3 (Figures 4C–E, Supplementary Figure S3). Consistently, the t-SNE confirmed that only cluster1, cluster3, and cluster4 could be significantly separated (Figure 4F). In total, 233 sepsis patients were classified into three subtypes, including cluster1 (n = 144), cluster 3 (n = 26), and cluster4 (n = 39) based on gene expression levels, which were chosen for subsequent analysis.

FIGURE 4

Identification of the DEGs Between Different Clusters in Sepsis

The DEGs among different molecular subtypes of sepsis patients were identified based on whole genome expression profiles. Eight hundred and sixty-seven (473 up-regulated and 394 down-regulated) DEGs were identified between cluster 1 and cluster 3 samples (Supplementary Figures S4A,B). The top 10 overexpressed and suppressed DEGs between cluster 1 and cluster 3 were exhibited in Table 2. In addition, a total of 420 (237 up-regulated and 183 down-regulated) DEGs were screened after comparing cluster 1 with cluster 4 samples (Supplementary Figures S4C,D). The top 10 overexpressed and suppressed DEGs between cluster 1 and cluster 4 were exhibited in Table 3. Moreover, one thousand four hundred and forty-five (696 up-regulated and 749 down-regulated) DEGs were acquired after comparing cluster 3 with cluster 4 samples (Supplementary Figures S4E,F). The top 10 overexpressed and suppressed DEGs between cluster 3 and cluster 4 were exhibited in Table 4. Overall, we identified 1803 specific genes among three subgroups.

GSVA Analysis in Different Molecular Subtypes of Sepsis

GSVA analysis was performed among distinct molecular subtypes of sepsis patients based on the expression profiles of DEGs. In total, 103 differentially enriched gene pathways, comprising 59 pathways with activation and 44 pathways with suppression were screened when comparing cluster1 with cluster3 samples (Supplementary Table S1). The top 20 differentially enriched gene pathways between cluster1 and cluster3 are exhibited in a heatmap (Figure 5A). A total of 90 differentially enriched gene pathways, consisting of 42 pathways with activation and 48 pathways with suppression, were obtained when comparing cluster1 with cluster4 samples (Supplementary Table S2). The top 20 differentially enriched gene pathways between cluster1 and cluster4 are exhibited in a heatmap (Figure 5B). A total of 110 differentially enriched gene pathways were obtained, including 56 pathways with activation and 54 pathways with suppression after comparing cluster3 with cluster4 samples (Supplementary Table S3). The top 20 differentially enriched gene pathways between cluster3 and cluster4 are exhibited in a heatmap (Figure 5C). Finally, a total of 34 enriched gene pathways among clusters were identified and presented in a heatmap (Figure 5D).

FIGURE 5

Identification and Enrichment Analysis of Co-DEGs

The DEGs identified by the DEGs and WGCNA methods in GSE154918 and GSE54514 datasets were intersected with the specific genes among different molecular subtypes of sepsis patients. Eventually, a total of 40 co-DEGs were screened (Figure 6A). Spearman’s correlation analysis was performed to illuminate the correlation patterns among these 40 co-DEGs, and the area of the pie chart represents the exact value of correlation coefficients. We found a significant correlation among these co-DEGs. For example, GPR84 show strong antagonistic effects with GYG1 (coefficient = 0.968), FCGR1B (coefficient = 0.968), and ANKRD22 (coefficient = 0.912). Simultaneously, CARD11 could also present synergistic effects with GYG1 (coefficient = -0.908), GPR84 (coefficient = -0.91), and BLOC1S1 (coefficient = -0.918) (Figure 6B). These results suggested the regulatory balance among these co-DEGs. In addition, the gene relationship network indicated that genes with a correlation coefficient of more than 0.9 were as follows: IRAK3, IL18R1, ATP6V1CA, GYG1, ATP9A, NMT2, ITK, EIF4B, CCND2, CARD11, FBXO21, BLOC1S1, GPR84, LRG1, ANKRD22, FCGR1B (Figure 6C), indicating that the molecular subtypes of sepsis may be the results of multi-gene interactions. Subsequently, we assessed the enrichment pathways in which these co-DEGs involve. The GO enrichment analysis suggested that these co-DEGs were primarily enriched in immune-related biological functions and pathways, such as immune response, T cell receptor signaling pathway, T cell costimulation, positive regulation of NF-kappaB transcription factor activity 4, positive regulation of T cell proliferation (Table 5 and Figures 6D,E). KEGG pathway enrichment analysis indicated that co-DEGs chiefly participated in autoimmune and inflammation-related diseases including Graft−versus−host disease, Inflammatory bowel disease, and Rheumatoid arthritis. Additionally, Immune responses and hematopoietic-related signaling pathways such as T cell receptor signaling pathway and Hematopoietic cell lineage were also closely associated with these co-DEGs (Figure 6F). Consistently, Reactome enrichment analysis also revealed that these co-DEGs were mainly implicated in immunomodulation and signal transduction such as Translocation of ZAP-70 to immunological synapse, Generation of second messenger molecules, Costimulation by the CD28 family, TCR signaling, and Downstream TCR signaling (Figure 6G). Combining these results, it can be inferred that abnormal immune function may be the critical pathogenesis of sepsis. In addition, the PPI network of 40 co-DEGs, including 22 nodes and 27 edges, was analyzed by the STRING website and visualized by the Cytoscape software (Supplementary Figure S5).

FIGURE 6

Identification of Potential Biomarkers Associated With Sepsis Using LASSO Regression

To further screen out potential biomarkers for sepsis, we performed the LASSO regression analysis based on the expression profile of co-DEGs. We separated a total of 303 samples in the combined dataset (GSE9960, GSE13904, and GSE54514 datasets, 233 sepsis and 70 control samples) into a training set (70%) and a validation set (30%). Eventually, a total of 25 potential biomarkers with a non-zero coefficient were obtained (Figures 7A,B). The AUC of the 25-gene signature was 0.9051 in the training set and 0.7955 in the validation set (Figure 7C), which suggests that the 25-gene-based model may be able to correctly diagnose sepsis. Subsequently, we also utilized external datasets to ascertain the diagnostic value of these potential biomarkers. The value of the AUC of these potential biomarkers in the GSE154918 and GSE69063 datasets were exhibited in a scatter plot (Figure 7D). Genes with AUC >0.95 in both the GSE154918 and GSE69063 datasets were as follows: ANKRD22, GPR84, GYG1, BLOC1S1, CARD11, NOG, and LRG1 (Figure 7E). These results reveal that the 7 core genes identified have the most ability to differentiate sepsis from healthy controls.

FIGURE 7

Verification of the Hub Genes Expression

To further demonstrate whether these 7 hub genes were worth using in clinical practice, we validated the expression of these core genes in different datasets. In GSE154918 and GSE69063 datasets, ANKRD22, GPR84, GYG1, BLOC1S1, and LRG1 were all highly expressed in sepsis patients, whereas NOG and CARD11 were dramatically decreased in sepsis patients (Figures 8A,B). Meanwhile, these 7 core genes were also subtype-specific biomarkers. Therefore, we speculated that these 7 hub genes have a high diagnostic ability in sepsis patients and are closely related to different subtypes of sepsis.

FIGURE 8

Discussion

Sepsis is a life-threatening inflammatory response syndrome caused by an unbalanced response of the host to various infection processes (van der Poll et al., 2017). Due to a lack of timely early diagnosis and treatment, sepsis has become one of the diseases with a high fatality and disability rate worldwide, and reducing the mortality rate has become the ultimate goal of its treatment. With in-depth studies of the pathogenesis of sepsis and the continuous attempt at treatment methods, distinct biomarkers have been employed in the diagnosis and treatment monitoring of sepsis (Faix, 2013; Nguyen et al., 2016). However, the interaction of the multiple genes involved in diverse biological functions may cause individual differences and complex pathophysiological mechanisms of sepsis (Maslove and Wong, 2014), which leads to unsatisfactory clinical diagnosis of a single biomarker. Therefore, it is urgent to identify distinct molecular subtypes of sepsis and elucidate the core genes and pathways associated with sepsis.

In our current study, we identified the DEGs and specific gene modules associated with sepsis in GSE154918 and GSE25504 datasets, utilizing the DEGs and WGCNA methods, respectively. Moreover, the “Consensus ClusterPlus” package was used to unsupervised cluster the patients of 233 sepsis patients from the GSE9960, GSE13904, and GSE54514 datasets. Afterward, a total of 40 co-DEGs were obtained by intersecting the DEGs, specific gene modules, and molecular cluster-related genes. GO, KEGG and Reactome enrichment analysis of these co-DEGs all indicated that sepsis was closely related to immune response and signal transduction, which may be the primary factor leading to sepsis progression. Therefore, screening of co-DEGs and identification of enrichment pathways significantly reduced the scope of our research, which may be able to ascertain more effective biomarkers for the early diagnosis and treatment of sepsis.

Defining different cut-off values plays a decisive role in DEG analysis. It is widely recognized that a more stringent screening criterion makes the results more convincing. However, according to the strict criterion for DEGs analysis (adjusted P-value < 0.05 and |logFC| > 1) (Liu et al., 2021; Lu et al., 2021; Su et al., 2021), we only obtained a total of 208 up-regulated and 380 down-regulated genes in the GSE154918 datasets, and 119 up-regulated and 302 down-regulated genes in the GSE25504 dataset. Therefore, we selected relatively rough cut-off values (p-value < 0.05 and |logFC| > 0.5) to identify as many DEGs and co-DEGs as possible. We finally acquired 3,294 DEGs, including 1,671 up-regulated and 1,623 down-regulated genes in the GSE154918 dataset, and 1738 DEGs, including 945 up-regulated and 793 down-regulated genes in the GSE25504 dataset.

Different sepsis patients have distinct prognoses, which may be caused by multiple pathways in which different key genes are involved. In our study, we performed an unsupervised cluster for sepsis patients based on whole-genome expression profiles and clustered 233 sepsis patients into three molecular subtypes. However, only 79 overlapped cluster-specific DEGs were screened when overlapping these DEGs among three subgroups. Therefore, to identify as many cluster-specific DEGs as possible, we merged these results and applied these 1803 merged cluster-specific DEGs for subsequent analysis. The consequences of GSVA analysis indicated that 103 differentially enriched gene pathways between cluster1 and 3 were screened, 90 differentially enriched gene pathways between clusters 1 and 4 were identified, and 110 differentially enriched gene pathways between clusters 1 and 4 were screened. These enriched gene pathways, such as metabolism, oxidative phosphorylation, autophagy regulation, and VEGF pathways have proved to be associated with the sepsis prognosis. Metabolism, the chemical reaction necessary for cell survival, usually retains the homeostasis between anabolism and catabolism under normal conditions (Qiu P. et al., 2019; Judge and Dodd, 2020; Xu et al., 2021). Additionally, metabolism involves the regulation of diverse cellular pathways, thus offering a large amount of energy for cells to ensure the execution of the function (Mulukutla et al., 2016; Panda, 2016). Oxidative phosphorylation is a coupling reaction that usually occurs in the mitochondrial inner membrane, which is characterized by the generation of cellular ATP based on mitochondrial electron transfer reactions (Ashton et al., 2018; Nolfi-Donegan et al., 2020; Głombik et al., 2021). Targeting mitochondrial oxidative phosphorylation is generally considered to be a novel strategy for the treatment of sepsis. Autophagy has been widely explored in numerous diseases including sepsis (Levy, 2007; Su et al., 2019). As previously described, autophagy exerts effects on immune regulation and inhibition of tissue damage after sepsis through regulating the expression of various immune cells (Oami et al., 2017). Moreover, autophagy has been demonstrated as an effective target for alleviating oxidative stress-induced organ failure after sepsis (Thiessen et al., 2017). VEGF pathways mainly participated in the regulation of blood vessel growth and play a vital role in promoting endothelial proliferation, migration, and survival by preserving the homeostasis of microvasculature (Ferrara and Gerber, 2001; Breen, 2007; Apte et al., 2019). With the progression of sepsis, VEGF activation results in vascular leak and dysfunction of host response, eventually leading to sepsis-related hypotension (van der Flier et al., 2005; Schuetz et al., 2011). To further screen sepsis-related hub genes, we intersected the DEGs identified by the DEGs and WGCNA methods with subtypes-related specific genes. A total of 40 co-DEGs were identified, most of which existed with significant correlation. The results of GO, KEGG, and Reactome enrichment analysis revealed that the immune response may be closely connected with the severity of sepsis. Therefore, we speculated that the interaction of immune response with autophagy, VEGF, oxidative stress, and metabolic pathways may be the major factor leading to the progression of sepsis (Minion and Tewari, 2018; Yin et al., 2019; McBride et al., 2020; Huff et al., 2021). These results, combined with our findings, indicate that these differentially activated pathways could develop potential therapeutic targets for sepsis patients with distinctive molecular subtypes.

In our current study, we received a 25-gene signature based on the LASSO model, which can accurately diagnose sepsis in both the train and validation datasets. Among them, ANKRD22, GPR84, GYG1, BLOC1S1, CARD11, NOG, and LRG1 were further identified as the most relevant key genes based on AUCs in the GSE154918 and GSE69063 datasets. Bioinformatics analysis demonstrated that ANKRD22 and NOG can serve as potential biomarkers for the progression of cancers (Tarragona et al., 2012; Qiu Y. et al., 2019; Wu et al., 2021). GPR84 is a kind of G-protein-coupled receptor activated by free fatty acids (FFA) and plays a critical role in regulating lipid metabolism (Paulsen et al., 2014; Recio et al., 2018). Enhanced GPR84 is closely related to the activation of inflammation, thus exacerbating the development of adiposity and diabesity (Nagasaki et al., 2012). GYG1 deficiency is known to be associated with polyglucosan body myopathy (Malfatti et al., 2014) and BLOC1S1 is widely recognized as a degradation substrate for IRE1alpha (Hur et al., 2012). Mutations in CARD1, a protein carrying special caspase-related recruitment domains, can result in the poor prognosis of diffuse large B-cell lymphoma (Dong et al., 2020). LRG1, which serves as a novel angiogenic factor, is required for the regulation of pathogenic angiogenesis (Wang et al., 2017). Although no previous studies have reported on the association between these 7 key genes and sepsis, the results of our analysis indicate that these hub genes may be potential markers for early diagnosis of sepsis.

To further explore and verify the clinical application value, we validated the expression levels of these 7 core genes in external datasets. All of the genes varied significantly between normal and sepsis samples in the GSE154918 and GSE69063 datasets. Additionally, sepsis with different molecular subtypes also exhibited the distinct expression of these 7 core genes. Further research needs to be carried out to elucidate the molecular mechanism of the 7 hub genes involved in the pathogenesis of sepsis molecular subtypes.

Several limitations have to be taken into account in our current study. First, the datasets used for analyses contained different sample sizes of normal and sepsis patients, which may influence the accuracy of the analytical results. Second, these datasets were downloaded from a publically available database, and lacked information on main clinical features such as sex, age, complications, recurrence rate, and individual therapeutic effect. A supplementary study is required, with a more detailed analysis of the demographic and clinical characteristics of sepsis. In addition, relatively rough cut-off values (p-value < 0.05 and |logFC| > 0.5) may influence the accuracy of the results. Moreover, the results of our analysis need to be confirmed in vitro, in vivo, and in clinical trials studies.

Conclusion

In conclusion, we identified 40 co-DEGs and several immune response pathways related to sepsis prognosis using various bioinformatics analyses. We constructed a 25-gene signature diagnostic model based on LASSO regression analysis, which has a high value for the early diagnosis of sepsis. There are remarkable differences in ANKRD22, GPR84, GYG1, BLOC1S1, CARD11, NOG, and LRG1 gene expression and enriched pathways among different molecular subgroups of sepsis, which may be the key factors leading to the heterogeneity of clinical symptoms and prognosis in patients with sepsis. Our current study provides novel diagnostic and therapeutic biomarkers for sepsis molecular subtypes.

References

• 1

  ApteR. S.ChenD. S.FerraraN. (2019). VEGF in signaling and disease: beyond discovery and development. Cell176, 1248–1264. 10.1016/j.cell.2019.01.021

  • CrossRef
  • Google Scholar

• 2

  AshtonT. M.McKennaW. G.Kunz-SchughartL. A.HigginsG. S. (2018). Oxidative phosphorylation as an emerging target in cancer therapy. Clin. Cancer Res.24, 2482–2490. 10.1158/1078-0432.CCR-17-3070

  • CrossRef
  • Google Scholar

• 3

  BaderG. D.HogueC. W. (2003). An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinforma.4, 2. 10.1186/1471-2105-4-2

  • CrossRef
  • Google Scholar

• 4

  BanwaitJ. K.BastolaD. R. (2015). Contribution of bioinformatics prediction in microRNA-based cancer therapeutics. Adv. Drug Deliv. Rev.81, 94–103. 10.1016/j.addr.2014.10.030

  • CrossRef
  • Google Scholar

• 5

  BarrettT.WilhiteS. E.LedouxP.EvangelistaC.KimI. F.TomashevskyM.et al (2013). NCBI GEO: archive for functional genomics data sets--update. Nucleic Acids Res.41, D991–D995. 10.1093/nar/gks1193

  • CrossRef
  • Google Scholar

• 6

  BreenE. C. (2007). VEGF in biological control. J. Cell. Biochem.102, 1358–1367. 10.1002/jcb.21579

  • CrossRef
  • Google Scholar

• 7

  ChengQ.ChenX.WuH.DuY. (2021). Three hematologic/immune system-specific expressed genes are considered as the potential biomarkers for the diagnosis of early rheumatoid arthritis through bioinformatics analysis. J. Transl. Med.19, 18. 10.1186/s12967-020-02689-y

  • CrossRef
  • Google Scholar

• 8

  DavenportE. E.BurnhamK. L.RadhakrishnanJ.HumburgP.HuttonP.MillsT. C.et al (2016). Genomic landscape of the individual host response and outcomes in sepsis: a prospective cohort study. Lancet. Respir. Med.4, 259–271. 10.1016/S2213-2600(16)00046-1

  • CrossRef
  • Google Scholar

• 9

  DongG.ChanudetE.ZengN.AppertA.ChenY. W.AuW. Y.et al (2020). A20, ABIN-1/2, and CARD11 mutations and their prognostic value in gastrointestinal diffuse large B-cell lymphoma. Clin. Cancer Res.17, 1440–1451. 10.1158/1078-0432.CCR-10-1859

  • CrossRef
  • Google Scholar

• 10

  FaixJ. D. (2013). Biomarkers of sepsis. Crit. Rev. Clin. Lab. Sci.50, 23–36. 10.3109/10408363.2013.764490

  • CrossRef
  • Google Scholar

• 11

  FerraraN.GerberH. P. (2001). The role of vascular endothelial growth factor in angiogenesis. Acta Haematol.106, 148–156. 10.1159/000046610

  • CrossRef
  • Google Scholar

• 12

  FleischmannC.ScheragA.AdhikariN. K.HartogC. S.TsaganosT.SchlattmannP.et al (2016). Assessment of global incidence and mortality of hospital-treated sepsis. Current estimates and limitations. Am. J. Respir. Crit. Care Med.193, 259–272. 10.1164/rccm.201504-0781OC

  • CrossRef
  • Google Scholar

• 13

  GłombikK.DetkaJ.BudziszewskaB. (2021). Hormonal regulation of oxidative phosphorylation in the brain in health and disease. Cells10, 2937. 10.3390/cells10112937

  • CrossRef
  • Google Scholar

• 14

  HuY.ChengX.QiuZ.ChenX. (2021). Identification of metabolism-associated molecular subtypes of chronic obstructive pulmonary disease. Int. J. Chron. Obstruct. Pulmon. Dis.16, 2351–2362. 10.2147/COPD.S316304

  • CrossRef
  • Google Scholar

• 15

  Huang daW.ShermanB. T.LempickiR. A. (2009). Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc.4, 44–57. 10.1038/nprot.2008.211

  • CrossRef
  • Google Scholar

• 16

  HuffL. A.YanS.ClemensM. G. (2021). Mechanisms of ataxia telangiectasia mutated (ATM) control in the DNA damage response to oxidative stress, epigenetic regulation, and persistent innate immune suppression following sepsis. Antioxidants (Basel)10, 1146. 10.3390/antiox10071146

  • CrossRef
  • Google Scholar

• 17

  HurK. Y.SoJ. S.RudaV.Frank-KamenetskyM.FitzgeraldK.KotelianskyV.et al (2012). IRE1α activation protects mice against acetaminophen-induced hepatotoxicity. J. Exp. Med.209, 307–318. 10.1084/jem.20111298

  • CrossRef
  • Google Scholar

• 18

  JudgeA.DoddM. S. (2020). Metabolism. Essays Biochem.64, 607–647. 10.1042/EBC20190041

  • CrossRef
  • Google Scholar

• 19

  LangfelderP.HorvathS. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinforma.9, 559. 10.1186/1471-2105-9-559

  • CrossRef
  • Google Scholar

• 20

  LevyR. J. (2007). Mitochondrial dysfunction, bioenergetic impairment, and metabolic down-regulation in sepsis. Shock28, 24–28. 10.1097/01.shk.0000235089.30550.2d

  • CrossRef
  • Google Scholar

• 21

  LiC.PanJ.LuoJ.ChenX. (2021). Prognostic characterization of immune molecular subtypes in non-small cell lung cancer to immunotherapy. BMC Pulm. Med.21, 389. 10.1186/s12890-021-01765-3

  • CrossRef
  • Google Scholar

• 22

  LiuT.LiX.CuiY.MengP.ZengG.WangY.et al (2021). Bioinformatics analysis identifies potential ferroptosis key genes in the pathogenesis of intracerebral hemorrhage. Front. Neurosci.15, 661663. 10.3389/fnins.2021.661663

  • CrossRef
  • Google Scholar

• 23

  LuT.XuR.LiQ.ZhaoJ. Y.PengB.ZhangH.et al (2021). Systematic profiling of ferroptosis gene signatures predicts prognostic factors in esophageal squamous cell carcinoma. Mol. Ther. Oncolytics21, 134–143. 10.1016/j.omto.2021.02.011

  • CrossRef
  • Google Scholar

• 24

  MalfattiE.NilssonJ.Hedberg-OldforsC.Hernandez-LainA.MichelF.Dominguez-GonzalezC.et al (2014). A new muscle glycogen storage disease associated with glycogenin-1 deficiency. Ann. Neurol.76, 891–898. 10.1002/ana.24284

  • CrossRef
  • Google Scholar

• 25

  MarkwartR.CondottaS. A.RequardtR. P.BorkenF.SchubertK.WeigelC.et al (2015). Immunosuppression after sepsis: systemic inflammation and sepsis induce a loss of nave T-cells but No enduring cell-autonomous defects in T-cell function. PLoS One9 (12), e115094. 10.1371/journal.pone.0115094

  • CrossRef
  • Google Scholar

• 26

  MasloveD. M.WongH. R. (2014). Gene expression profiling in sepsis: timing, tissue, and translational considerations. Trends Mol. Med.20, 204–213. 10.1016/j.molmed.2014.01.006

  • CrossRef
  • Google Scholar

• 27

  McBrideM. A.OwenA. M.StothersC. L.HernandezA.LuanL.BurelbachK. R.et al (2020). The metabolic basis of immune dysfunction following sepsis and trauma. Front. Immunol.11, 1043. 10.3389/fimmu.2020.01043

  • CrossRef
  • Google Scholar

• 28

  MinionL. E.TewariK. S. (2018). Cervical cancer - state of the science: from angiogenesis blockade to checkpoint inhibition. Gynecol. Oncol.148, 609–621. 10.1016/j.ygyno.2018.01.009

  • CrossRef
  • Google Scholar

• 29

  MulukutlaB. C.YongkyA.LeT.MashekD. G.HuW. S. (2016). Regulation of glucose metabolism - a perspective from cell bioprocessing. Trends Biotechnol.34, 638–651. 10.1016/j.tibtech.2016.04.012

  • CrossRef
  • Google Scholar

• 30

  NagasakiH.KondoT.FuchigamiM.HashimotoH.SugimuraY.OzakiN.et al (2012). Inflammatory changes in adipose tissue enhance expression of GPR84, a medium-chain fatty acid receptor: TNFα enhances GPR84 expression in adipocytes. FEBS Lett.586, 368–372. 10.1016/j.febslet.2012.01.001

  • CrossRef
  • Google Scholar

• 31

  NguyenH. B.JaehneA. K.JayaprakashN.SemlerM. W.HegabS.YatacoA. C.et al (2016). Early goal-directed therapy in severe sepsis and septic shock: insights and comparisons to ProCESS, ProMISe, and ARISE. Crit. Care20, 160. 10.1186/s13054-016-1288-3

  • CrossRef
  • Google Scholar

• 32

  Nolfi-DoneganD.BraganzaA.ShivaS. (2020). Mitochondrial electron transport chain: oxidative phosphorylation, oxidant production, and methods of measurement. Redox Biol.37, 101674. 10.1016/j.redox.2020.101674

  • CrossRef
  • Google Scholar

• 33

  NovosadS. A.SapianoM. R.GriggC.LakeJ.RobynM.DumyatiG.et al (2016). Vital signs: epidemiology of sepsis: prevalence of health care factors and opportunities for prevention. MMWR. Morb. Mortal. Wkly. Rep.65, 864–869. 10.15585/mmwr.mm6533e1

  • CrossRef
  • Google Scholar

• 34

  OamiT.WatanabeE.HatanoM.SunaharaS.FujimuraL.SakamotoA.et al (2017). Suppression of T Cell autophagy results in decreased viability and function of T cells through accelerated apoptosis in a murine sepsis model. Crit. Care Med.45, e77–e85. 10.1097/CCM.0000000000002016

  • CrossRef
  • Google Scholar

• 35

  PandaS. (2016). Circadian physiology of metabolism. Science354, 1008–1015. 10.1126/science.aah4967

  • CrossRef
  • Google Scholar

• 36

  PaulsenS. J.LarsenL. K.HansenG.ChelurS.LarsenP. J.VrangN.et al (2014). Expression of the fatty acid receptor GPR120 in the gut of diet-induced-obese rats and its role in GLP-1 secretion. PLoS One9, e88227. 10.1371/journal.pone.0088227

  • CrossRef
  • Google Scholar

• 37

  PierrakosC.VincentJ. L. (2010). Sepsis biomarkers: a review. Crit. Care14, R15. 10.1186/cc8872

  • CrossRef
  • Google Scholar

• 38

  QiuP.LiuY.ZhangJ. (2019). Review: the role and mechanisms of macrophage autophagy in sepsis. Inflammation42, 6–19. 10.1007/s10753-018-0890-8

  • CrossRef
  • Google Scholar

• 39

  QiuY.YangS.PanT.YuL.LiuJ.ZhuY.et al (2019). ANKRD22 is involved in the progression of prostate cancer. Oncol. Lett.18, 4106–4113. 10.3892/ol.2019.10738

  • CrossRef
  • Google Scholar

• 40

  RaoJ.LiW.ChenC. (2021). Pyroptosis-mediated molecular subtypes and tumor microenvironment infiltration characterization in colon cancer. Front. Cell Dev. Biol.9, 766503. 10.3389/fcell.2021.766503

  • CrossRef
  • Google Scholar

• 41

  RecioC.LucyD.PurvisG. S. D.IvesonP.ZeboudjL.IqbalA. J.et al (2018). Activation of the immune-metabolic receptor GPR84 enhances inflammation and phagocytosis in macrophages. Front. Immunol.9, 1419. 10.3389/fimmu.2018.01419

  • CrossRef
  • Google Scholar

• 42

  RelloJ.Valenzuela-SánchezF.Ruiz-RodriguezM.MoyanoS. (2017). Sepsis: a review of advances in management. Adv. Ther.34, 2393–2411. 10.1007/s12325-017-0622-8

  • CrossRef
  • Google Scholar

• 43

  RitchieM. E.PhipsonB.WuD.HuY.LawC. W.ShiW.et al (2015). Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res.43, e47. 10.1093/nar/gkv007

  • CrossRef
  • Google Scholar

• 44

  SandquistM.WongH. R. (2014). Biomarkers of sepsis and their potential value in diagnosis, prognosis and treatment. Expert Rev. Clin. Immunol.10, 1349–1356. 10.1586/1744666X.2014.949675

  • CrossRef
  • Google Scholar

• 45

  SchuetzP.JonesA. E.AirdW. C.ShapiroN. I. (2011). Endothelial cell activation in emergency department patients with sepsis-related and non-sepsis-related hypotension. Shock36, 104–108. 10.1097/SHK.0b013e31821e4e04

  • CrossRef
  • Google Scholar

• 46

  SeymourC. W.KennedyJ. N.WangS.ChangC. H.ElliottC. F.XuZ.et al (2019). Derivation, validation, and potential treatment implications of novel clinical phenotypes for sepsis. JAMA321, 2003–2017. 10.1001/jama.2019.5791

  • CrossRef
  • Google Scholar

• 47

  ShiX.CaoY.ZhangX.GuC.LiangF.XueJ.et al (2021). Comprehensive analysis of N6-methyladenosine RNA methylation regulators expression identify distinct molecular subtypes of myocardial infarction. Front. Cell Dev. Biol.9, 756483. 10.3389/fcell.2021.756483

  • CrossRef
  • Google Scholar

• 48

  SuL. J.ZhangJ. H.GomezH.MuruganR.HongX.XuD.et al (2019). Reactive oxygen species-induced lipid peroxidation in apoptosis, autophagy, and ferroptosis. Oxid. Med. Cell. Longev.2019, 5080843. 10.1155/2019/5080843

  • CrossRef
  • Google Scholar

• 49

  SuW.ZhaoY.WeiY.ZhangX.JiJ.YangS.et al (2021). Exploring the pathogenesis of psoriasis complicated with atherosclerosis via microarray data analysis. Front. Immunol.12, 667690. 10.3389/fimmu.2021.667690

  • CrossRef
  • Google Scholar

• 50

  SzklarczykD.FranceschiniA.WyderS.ForslundK.HellerD.Huerta-CepasJ.et al (2015). STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res.43, D447–D452. 10.1093/nar/gku1003

  • CrossRef
  • Google Scholar

• 51

  TarragonaM.PavlovicM.Arnal-EstapéA.UrosevicJ.MoralesM.GuiuM.et al (2012). Identification of NOG as a specific breast cancer bone metastasis-supporting gene. J. Biol. Chem.287, 21346–21355. 10.1074/jbc.M112.355834

  • CrossRef
  • Google Scholar

• 52

  ThiessenS. E.DereseI.DerdeS.DufourT.PauwelsL.BekhuisY.et al (2017). The role of autophagy in critical illness-induced liver damage. Sci. Rep.7, 14150. 10.1038/s41598-017-14405-w

  • CrossRef
  • Google Scholar

• 53

  van der FlierM.van LeeuwenH. J.van KesselK. P.KimpenJ. L.HoepelmanA. I.GeelenS. P.et al (2005). Plasma vascular endothelial growth factor in severe sepsis. Shock23, 35–38. 10.1097/01.shk.0000150728.91155.41

  • CrossRef
  • Google Scholar

• 54

  van der PollT.van de VeerdonkF. L.SciclunaB. P.NeteaM. G. (2017). The immunopathology of sepsis and potential therapeutic targets. Nat. Rev. Immunol.17, 407–420. 10.1038/nri.2017.36

  • CrossRef
  • Google Scholar

• 55

  WangY.XuJ.ZhangX.WangC.HuangY.DaiK.et al (2017). TNF-α-induced LRG1 promotes angiogenesis and mesenchymal stem cell migration in the subchondral bone during osteoarthritis. Cell Death Dis.8, e2715. 10.1038/cddis.2017.129

  • CrossRef
  • Google Scholar

• 56

  WilkersonM. D.HayesD. N. (2010). ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics26, 1572–1573. 10.1093/bioinformatics/btq170

  • CrossRef
  • Google Scholar

• 57

  WuY.LiuH.GongY.ZhangB.ChenW. (2021). ANKRD22 enhances breast cancer cell malignancy by activating the Wnt/β-catenin pathway via modulating NuSAP1 expression. Bosn. J. Basic Med. Sci.21, 294–304. 10.17305/bjbms.2020.4701

  • CrossRef
  • Google Scholar

• 58

  XuW.JanochaA. J.ErzurumS. C. (2021). Metabolism in pulmonary hypertension. Annu. Rev. Physiol.83, 551–576. 10.1146/annurev-physiol-031620-123956

  • CrossRef
  • Google Scholar

• 59

  YinX.XinH.MaoS.WuG.GuoL. (2019). The role of autophagy in sepsis: protection and injury to organs. Front. Physiol.10, 1071. 10.3389/fphys.2019.01071

  • CrossRef
  • Google Scholar

• 60

  YuG.WangL. G.HanY.HeQ. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS16, 284–287. 10.1089/omi.2011.0118

  • CrossRef
  • Google Scholar

• 61

  ZhangS.LiuF.WuZ.XieJ.YangY.QiuH. (2020). Contribution of m6A subtype classification on heterogeneity of sepsis. Ann. Transl. Med.8, 306. 10.21037/atm.2020.03.07

  • CrossRef
  • Google Scholar