Identification of an Immune-Related Gene Diagnostic Model and Potential Drugs in Sepsis Using Bioinformatics and Pharmacogenomics Approaches

Introduction

Sepsis is described as the syndrome consisting of complex biochemical and pathophysiological dysregulation.¹ Post-sepsis syndrome will increase the risk for infections and the incidence of mental health problems, cognitive impairment, cardiovascular events, and renal failure.^2–5 Apart from antimicrobial agents, some alternative approach therapies have shown some additive benefits on survival in patients with sepsis, such as host defense peptides that not only directly kill pathogens, but also possess immunomodulatory properties.⁶ Under the combined treatment of one or more therapies, although the mortality rate of sepsis patients has significantly decreased, it is still as high as 20%. In addition, early diagnosis is difficult to achieve in everyday clinical practice. Therefore, understanding the pathological mechanism and identifying potential diagnostic and therapeutic markers is pivotal to improve clinical outcomes and reduce mortality of sepsis.^7,8

Generally, the pathophysiology of sepsis is regarded as an initial hyper-inflammatory phase followed by the immunosuppressive phase.⁹ After inflammation, the immune system leads to the hypo-inflammatory state. Besides, the mass depletion of immune cells will leave patients vulnerable to secondary infection.^10–15 Numerous studies have demonstrated that the main factor of sepsis is the host’s response to infection.¹⁶ Hence, it is necessary to understand the role of immune cells in sepsis.

In this study, some public microarray datasets were comprehensively analyzed to evaluate the immune score of sepsis and identify differentially expressed genes (DEGs). Those hub DEGs related to immune score of sepsis were identified through weighted gene co-expression network analysis (WGCNA) and protein–protein interaction (PPI) network. Based on these immune-related hub genes, the diagnostic model for predicting sepsis was constructed based on 3 classification models, including Random Forest (RF), Support Vector Machine (SVM) and Decision Tree (DT).

Methods

Dataset Collection

All data were downloaded from the Gene Expression Omnibus (GEO) dataset. Datasets were searched with the keyword “Sepsis” and “homo sapiens”. Finally, a total of 3 datasets were collected, including GSE69063, GSE154918 and GSE28750. In these datasets, peripheral blood samples from adult patients were collected for analysis. The detailed information of these datasets is listed in Table 1. The integrated data of GSE69063 and GSE154918 datasets was used for the training set. GSE28750 dataset was used for the validation set. In GSE69063 and GSE154918 datasets, gene expression lists were downloaded to obtain the intersection genes. The ComBat function of SVA package was used to remove batch effects.

Table 1 Detail Information of 3 Datasets in Sepsis

Immunoinfiltration Analysis and Identification of DEGs in Sepsis

For the integrated data of GSE69063 and GSE154918 datasets, the ESTIMATE software was used to evaluate stromal, immune, and ESTIMATE scores. Single-sample gene set enrichment analysis (ssGSEA) was used to assess the score of 23 immune cells in sepsis. After data preprocessing, DEGs were screened using limma packages under the threshold value of false discovery rate (FDR) <0.05 and |Fold change (FC)| >1. Bidirectional hierarchical cluster analysis was performed using gene expression values based on Euclidean distance. The results were presented in the form of heat maps. For the training set, the rank sum test was utilized to analyze the expression of DEGs in sepsis. The expression of DEGs was also verified in the validation set. The box plot was used to show the gene expression and the relationship between DEGs and differential immune cell was analyzed.

WGCNA in Sepsis

WGCNA is a typical systematic bio-algorithm for describing the correlation profiles between gene expressions to construct gene co-expression networks. Cluster analysis was conducted for the modules. The Pearson method was utilized to analyze the correlation with sepsis. The module with the highest correlation with immune score was chosen as the hub module.

Identification of Immune-Related Key Genes and Protein–Protein Interaction (PPI) Analysis in Sepsis

The immune-related genes were obtained from import database. Immune-related key genes were identified by intersection with DEGs and genes in hub module of WGCNA. Additionally, PPI networks were constructed for the immune-related key genes using string database. Four algorithms, Maximal Clique Centrality (MCC), Maximum Neighborhood Component (MNC), degree and closeness, were used to identify hub genes. Intersection genes of key genes were identified as the final hub genes.

Construction of Diagnostic Model Based on Immune-Related Hub Genes

GSE69063 and GSE154918 datasets were used as the training set, and the GSE28750 dataset was used as the validation set. Immune-related hub genes were used as the features in the training set. The corresponding expression profiles of these immune-related hub genes were obtained. Three classification models were constructed, including RF, SVM and DT. The area under the receiver operating characteristic (ROC) curve (AUC) of the model was analyzed. The receiver operating characteristic is a comprehensive index reflecting the false-positive rate and true positive rate of continuous variables and the area under the curve (AUC) is an evaluation index of model performance.

Functional Analysis and Drug Prediction of Immune-Related Hub Genes in Sepsis

To further explore the biological functions of immune-related hub genes, David software was applied for functional analysis, including gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). Moreover, hub genes were selected to screen out drugs related to their proteins.

Expression Validation of Immune-Related Hub Genes in Sepsis by Reverse Transcriptase Polymerase Chain Reaction (RT-PCR)

From April 10 to July 16, 2022, the sepsis patients and normal individual in the General Hospital of Northern Theater Command PLA were recruited to perform RT-PCR to validate the expression of immune-related hub genes in blood samples from sepsis patients. The normal individuals were from the population undergoing routine physical examination at our hospital. The inclusion criteria of sepsis patients were as follows: (1) patients were diagnosed with sepsis based on the criteria of the American Society of Intensive Care Medicine (SCCM), the European Society of Intensive Care Medicine (ESICM), the American College of Chest Physicians (ACCP), and International sepsis;^1,17 (2) patients had a possible source of infection and any two of the following indicators including body temperature above 38°C or below 36°C, heart rate exceeding 90 beats per minute or 2 standard deviations above the normal heart rate range for different ages, respiration rate greater than 20 beats per minute or PaCO₂ less than 32 mmHg, white blood cell counts greater than 12.0 ×10⁹/L or less than 4.0 ×10⁹/L, or more than 10% of immature neutrophils; (3) patients had complete clinical data. The exclusion criteria of sepsis patients were as follows: (1) patients with septic shock, a history of cancer or other diseases (such as severe trauma, shock, and severe acute pancreatitis), chemotherapy, radiation, and so on. (2) patients had incomplete clinical data. The inclusion criteria of healthy individuals were as follows: (1) these individuals matched by gender and age; (2) no disease occurred in the two weeks before and after sampling. The exclusion criteria of healthy individuals were as follows: (1) those individuals who took glucocorticoids two weeks prior to sampling; (2) those individuals with a history of febrile illness, or any chronic/acute illness mildly associated with inflammation in the two weeks prior to sampling. Based on the above inclusion and exclusion criteria, a total of 8 sepsis patients and 7 healthy individuals were enrolled in the present study. Clinical information of 8 sepsis patients and 7 healthy individuals are listed in Table 2, with age and gender matching between the two groups (Supplementary Table 1). The blood samples of these individuals were collected for RT-PCR. GAPDH and ACTB were used as internal references.

Table 2 Clinical Information of 8 Sepsis Patients and 7 Healthy Individuals in RT-PCR

Statistical Analysis

T test was used for statistical analysis in the expression validation in RT-PCR. P < 0.05 was considered as statistical difference.

Results

Data Preprocessing

The gene expression list was downloaded from GSE69063 and GSE154918 datasets. A total of common 17,564 genes were obtained in two datasets. The principle component analysis (PCA) diagrams and box plots after removing batch effect are shown in Figure 1.

Figure 1 Box plots of gene data distribution in sepsis. (A) PCA diagrams before removing batch effect; (B) PCA diagrams after removing batch effect; (C) box plots before removing batch effect; (D) box plots after removing batch effect.

Immunoinfiltration Analysis in Sepsis

Immunoinfiltration analysis was performed on the training set. The score of 23 immune cells was evaluated using ssGSEA (Figure 2A). Compared to normal control, the infiltration of activated dendritic cell, CD56bright natural killer cell, gamma delta T cell, macrophage, mast cell, monocyte, neutrophil, regulatory T cell and type 17 T helper cell were significantly increased in sepsis. The infiltration of activated B cell, activated CD8 T cell, CD56dim natural killer cell, immature B cell, immature dendritic cell, natural killer cell, natural killer T cell, T follicular helper cell, and Type 1 T helper cell were remarkably reduced in sepsis. In addition, stromal score, immune score, and ESTIMATE score (Figure 2B-D) were evaluated. The stromal score and ESTIMATE score were higher, while immune score was lower in sepsis compared to normal controls.

Figure 2 Immunoinfiltration analysis in sepsis. (A) immune score of 23 immune cells. *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001. (B) stromal score; (C) immune score; (D) ESTIMATE score. Abbreviation: ns, not significant.

WGCNA in Sepsis

All genes from GSE69063 and GSE154918 were selected to construct a co-expression network. Due to sensitivity of WGCNA to abnormal samples, strict quality control procedures were performed to ensure the highest quality levels. The soft threshold was set to 4 (Figure 3A), a total of 13 modules were identified (Figure 3B). Significantly, turquoise module had the highest correlation with immune score and the most significant correlation with disease (Figure 3C). Therefore, the turquoise module (involving 12,459 genes) was considered as the hub module.

Figure 3 WGCNA in sepsis. (A) scale-free fitting index of different soft threshold power and mean connectivity of various soft threshold powers; (B) genes are divided into different modules; (C) heat map of the correlation between modular characteristic genes and sepsis.

Screening of Immune-Related Key Genes in Sepsis

Totally, 561 DEGs and 1793 immune-related genes were identified in sepsis. The volcano plot and heat map of top 100 DEGs are presented in Figure 4A and B, respectively. A total of 93 immune-related key genes were identified between 561 DEGs, 1793 immune-related genes and 12,459 genes in the hub module (Figure 4C).

Figure 4 Screening of immune-related key genes in sepsis. (A) volcano plot of top 100 DEGs; (B) heat map of top 100 DEGs; (C) Venn diagram of DEGs, genes in hub module, and immune-related genes.

PPI Analysis of Immune-Related Key Genes in Sepsis

The PPI network was constructing based on 93 immune-related key genes (Figure 5A). After screening through 4 algorithms (MCC, MNC, degree and closeness), 5 hub genes were identified (Figure 5B), including interleukin 7 receptor (IL7R), interleukin 10 (IL10), CD40 ligand (CD40LG), CD28 molecule (CD28) and lipocalin 2 (LCN2). The expression levels of above 5 genes were tested in training set (Figure 6A) and the validation set (Figure 6B). IL10 and LCN2 were significantly up-regulated, while IL7R, CD40LG and CD28 were remarkably down-regulated in sepsis. Furthermore, the association between these genes and immune cell was analyzed (Figure 6C). LCN2 was significantly positively related to activated dendritic cell, IL7R and CD28 were significantly negatively related to activated dendritic cell, and IL10 was significantly negatively correlated with immature B cell.

Figure 5 PPI analysis of immune-related key genes in sepsis. (A) PPI network of 93 immune-related key genes; (B) Venn diagram of hub genes.

Figure 6 The expression and relationship with immune cell analysis of 5 hub immune-related genes. (A) gene expression in the training set, ****p < 0.0001; (B) gene expression in the validation set, ****p < 0.0001; (C) relationship between genes and immune cell, *p < 0.05; **p < 0.01; ***p < 0.001.

Construction and Validation of Diagnostic Model Based on Immune-Related Hub Genes

Based on immune-related hub genes, three models (RF, SVM and DT) were established in training set (Figure 7A) and the validation set (Figure 7B). The AUC of the three models was high, and the AUC of SVM was the highest. In addition, the ROC of single IL7R, IL10, CD40LG, CD28 and LCN2 was analyzed in training set (Figure 8A) and the validation set (Figure 8B). The AUC of these genes was all more than 0.8. The diagnostic prediction model could effectively distinguish sepsis samples from normal samples. The differential expression analysis of immune-related hub genes between survivor vs non-survivor and shock vs non-shock in sepsis patients showed that LCN2 showed significant differences between shock and non-shock, while other genes did not show significant differences between the two groups (Figure 9).

Figure 7 The ROC curves of the three models in training set (A) and validation set (B).

Figure 8 The ROC curves of IL7R, IL10, CD40LG, CD28 and LCN2 in training set (A) and validation set (B).

Figure 9 The gene expression of immune-related genes between sepsis patients. (A) the gene expression between control, sepsis survivor and non-survivor patients; (B) the gene expression between control, shock and non-shock in sepsis patients; (C) the gene expression between sepsis survivor and non-survivor patients; (D) the gene expression between sepsis shock and non-shock patients; **p < 0.01; ****p < 0.0001. ns: not significant.

Enrichment Analysis and Drug Prediction of Immune-Related Hub Genes in Sepsis

According to GO analysis of 93 immune-related genes, immune response, external side of plasma membrane and C–C chemokine receptor activity was, respectively, the most significantly enriched biological process (Figure 10A), cytological component (Figure 10B) and molecular function (Figure 10C). Some significantly enriched signaling pathways were identified (Figure 10D), such as T cell receptor signaling pathway (involving CD28), NF-kappa B signaling pathway (involving CD40LG) and JAK-STAT signaling pathway (involving IL10). In addition, the DGIdb database was applied to find drugs related to proteins of IL7R, IL10, CD40LG, CD28 and LCN2. The results were imported into Cytoscape to build the relationship network (Figure 11), of which 17 candidate drugs might interact with IL 10, including CYCLOSPORINE. Six candidate drugs might interact with CD28 and 11 with CD40LG, and it is worth noting that CD40LG and CD28 were drug targets of ALDESLEUKIN.

Figure 10 Enrichment analyses of 93 immune-related genes in sepsis. (A) biological process (BP); (B) cytological component (CC); (C) molecular function (MF); (D) signaling pathway.

Figure 11 Regulatory networks between drugs and proteins encoded by IL7R, IL10, CD40LG, CD28 and LCN2. Red and blue color represents protein and drug, respectively.

Expression Validation of Immune-Related Hub Genes in Sepsis

Totally, 5 immune-related hub genes (IL10, LCN2, IL7R, CD40LG, and CD28) were used for expression validation by RT-PCR (Figure 12). IL10 and LCN2 were up-regulated, while IL7R, CD40LG and CD28 were down-regulated in the blood samples of sepsis patients. The expression trend consisted of the informatics analysis results.

Figure 12 Expression validations of immune-related hub genes in sepsis.

Discussion

In this study, a total of 5 hub genes were identified, including up-regulated IL10 and LCN2, and down-regulated IL7R, CD40LG and CD28. Significantly, diagnostic prediction model based on above 5 genes could be considered as reliable biomarkers for sepsis diagnosis.

IL10, associated with Gram-negative infections, has a negative effect in controlling bacterial growth and dissemination.¹⁸ Dreschers S et al found that IL10 activated CBMФ failed to increase glycolysis and showed reduced activation of the mechanistic target of rapamycin kinase (mTOR) pathway, which is important for survival in sepsis.¹⁹ A meta-analysis showed that IL10 gene polymorphisms (rs1800871, rs1800872, rs1800896) were associated with sepsis.²⁰ Studies have shown that IL-10 as a temporal regulator of the transition from early reversible sepsis to the late phase of irreversible shock.²¹ Additionally, NK cell-derived IL 10, produced in response to IL 15, is relevant to clinical manifestations in septic patients and critical for survival during sepsis.²² LCN2, originally isolated from neutrophil granules, is functionally associated with inflammation.^23,24 LCN2 was up-regulated in myeloid-derived suppressor cells from patients with severe sepsis or septic shock,²⁵ which was consistent with our study (Figure 9). In addition, LCN2 is related to organ failure degree and mortality in sepsis, and has a diagnostic value in sepsis-related acute respiratory distress syndrome.^26–28 Up-regulation of IL10 and LCN2 may play important roles in the inflammatory reaction of sepsis.

The mRNA level of IL7R is significantly reduced in sepsis.²⁹ In T cells, the expression increased of IL7R serves to enhance the effect of IL7 and improves survival in sepsis.³⁰ CD40 ligand transfers signals through CD40 for B-cell activation.^31,32 CD40LG is down-regulated in community acquired pneumonia patients presenting with organ failure.³³ CD28 plays important roles in regulating the activation of T lymphocytes, which are critical to combat infections.^34–36 A significantly decreased expression of CD28 was found in severe sepsis patients.^37,38 It is suggested that IL7R, CD40LG and CD28 may be associated with sepsis infection.

According to immunoinfiltration analysis, the infiltration of activated dendritic cell and immature B cell was significantly increased and reduced in sepsis. As a central link between adaptive and innate immune response, CD11c positive dendritic cells are critically associated with sepsis pathogenesis.³⁹ It is reported that myeloid cell-derived prostaglandins contribute to infection-associated apoptosis of immature B cells.⁴⁰ Notably, LCN2 was significantly positively related to activated dendritic cell, IL7R and CD28 were significantly negatively related to activated dendritic cell, and IL10 was significantly negatively correlated with immature B cell. LCN2 plays a significant function in the immune response, limiting bacterial growth and regulating innate immune cells.^41,42 IL7R could develop into plasma like dendritic cells.⁴³ The ligand expression of CD28 is necessary for the completion of cognate T-cell activation and the antigen presentation process. Thus, it can be seen that the interactive relationship of LCN2/IL7R/CD8-activated dendritic cell and IL10-immature B cell play crucial roles in the immune response of sepsis.

Based on KEGG analysis of 93 immune-related genes, four significantly enriched signaling pathways were identified, including T cell receptor signaling pathway (involving CD28), NF-kappa B signaling pathway (involving CD40LG) and JAK-STAT signaling pathway (involving IL10). Previous studies have shown that these signaling pathways participate in the physiological and pathological mechanisms of sepsis by influencing the immune response. For example, it has been demonstrated that the family of T cell co-inhibitory receptors inhibit CD8⁺ T cell functionality and contribute to post-septic immune incompetence.⁴⁴ The NF-kappa B family is involved in functions of inflammation, immunity, and stress response. In septic brains, inhibition of NF-kappa B signaling pathway enhances the completion of autophagy with a neuro-protective function.⁴⁵ It is found that expression and methylation alterations associated with JAK2-STAT pathway are relevant to sepsis patients infected with gram-negative bacteria.⁴⁶ Therefore, the key genes involved in these signaling pathways may be potential drug targets for sepsis treatment.

In addition, the DGIdb database was applied to find drugs related to proteins of IL7R, IL10, CD40LG, CD28 and LCN2. IL10 was a drug target of CYCLOSPORINE. CD40LG and CD28 were drug targets of ALDESLEUKIN. Cyclosporine A (CsA) has achieved some application effects in experimental acute inflammation such as sepsis, traumatic/hemorrhagic shock, and ischemia/reperfusion injury.⁴⁷ J Larche et al found that treatment with CYCLOSPORINE analogs could prevent sepsis-induced organ damage, decrease sepsis-related contractile dysfunction, and reduce the mortality rate through inhibition of caspase-3-like activity.⁴⁸ ALDESLEUKIN is an immunosuppressive drug that not only mobilizes immune effector cells with tumor killing ability but also alters immune response.^49–51 Research has shown that low-dose ALDESLEUKIN administration can induce selective amplification of regulatory T cells (Tregs) in mice and humans, thereby preventing the onset of autoimmune and inflammatory diseases.⁵² Our study indicated that CYCLOSPORINE and ALDESLEUKIN could be, respectively, used in clinical therapy of sepsis by targeting IL10, CD40LG and CD28.

Our study also has some limitations. Firstly, the sample size is relatively small, larger numbers of samples are further needed for expression validation of identified genes and potential drugs. Secondly, targeted studies on identified 5 diagnose/immune-related genes are needed to detect the expression of these genes in patients with sepsis. Thirdly, the incidence, physiology, and genetics of sepsis are known to differ between the sexes,^53,54 which require gender stratified sensitivity analysis. Finally, the relation between these gene expressions and the sepsis severity needs to be investigated in animal model.

Conclusion

In conclusion, our study identified a 5-gene diagnostic model based on immune infiltration analysis in sepsis, including IL7R, IL10, CD40LG, CD28 and LCN2. Two interactive relationship pairs were identified based on above genes, including interactive relationship of LCN2/IL7R/CD8-activated dendritic cell and IL10-immature B cell. Additionally, IL10 was a drug target of CYCLOSPORINE. CD40LG and CD28 were drug targets of ALDESLEUKIN. Our study may provide a new field in understanding the immunological mechanism of sepsis.