Gene Expression Scoring of Immune Activity Levels for Precision Use of Hydrocortisone in Vasodilatory Shock.

PURPOSE: Among patients with vasodilatory shock, gene expression scores may identify different immune states. We aimed to test whether such scores are robust in identifying patients' immune state and predicting response to hydrocortisone treatment in vasodilatory shock.
MATERIALS AND METHODS: We selected genes to generate continuous scores to define previously established subclasses of sepsis. We used these scores to identify a patient's immune state. We evaluated the potential for these states to assess the differential effect of hydrocortisone in two randomized clinical trials of hydrocortisone versus placebo in vasodilatory shock.
RESULTS: We initially identified genes associated with immune-adaptive, immune-innate, immune-coagulant functions. From these genes, 15 were most relevant to generate expression scores related to each of the functions. These scores were used to identify patients as immune-adaptive prevalent (IA-P) and immune-innate prevalent (IN-P). In IA-P patients, hydrocortisone therapy increased 28-day mortality in both trials (43.3% vs 14.7%, P = 0.028) and (57.1% vs 0.0%, P = 0.99). In IN-P patients, this effect was numerically reversed.
CONCLUSIONS: Gene expression scores identified the immune state of vasodilatory shock patients, one of which (IA-P) identified those who may be harmed by hydrocortisone. Gene expression scores may help advance the field of personalized medicine.

BACKGROUND

Vasodilatory shock is the most common form of shock in Intensive Care Unit (ICU) patients (1–3). Sepsis and severe burns are two common causes of vasodilatory shock and share similar pathophysiological mechanisms (3, 4). In their initial stages, these conditions are characterized by a pronounced systemic inflammatory response syndrome. Later, if immune homeostasis is not restored, both may develop a state of immunosuppression (3). Moreover, this immune response process is a complex interplay between innate and adaptive immune systems with an impact on the coagulation cascade (5, 6).

Given the presence of systemic inflammation (3), corticosteroids have been proposed as a therapeutic option in patients with vasodilatory shock. This is due to their immunomodulatory effects on both the adaptive and innate immune system and their successful application in diverse sets of inflammatory and auto-immune diseases (4, 5). However, the potential mortality benefits of corticosteroid therapy remain controversial due to conflicting clinical findings (6, 7). A plausible explanation for these varying results is the inherent clinical and biological heterogeneity of the immune response during sepsis and/or severe inflammation (5, 7, 8). For example, septic shock patients showing immunocompetent features may be harmed by hydrocortisone therapy due to its immunosuppressive effect on the adaptive immune system (9). Thus, identifying patient in different immune states may help target corticosteroid therapy.

Recently, several authors have attempted to identify more biologically homogeneous groups of sepsis patients (hereafter referred as sepsis subclasses) by applying various unsupervised learning techniques to both clinical and transcriptomic data (10–13). These subclasses may help separate underlying biological mechanisms and clinical characteristics (14).

We considered that previously established subclasses of sepsis could support the selection of genes to generate continuous scores that accurately reflect immune activity levels. These scores could then be used to stratify vasodilatory shock patients into those with adaptive or innate immune prevalent states. We hypothesized that patients with an immune-adaptive prevalent state would respond differently to hydrocortisone therapy compared with patients with an immune-innate prevalent state.

METHODS

Data source and participants

To identify unique biological subgroups of patients, we performed a retrospective study using a de-identified dataset pooling publicly available data from nine Gene Expression Omnibus (GEO) datasets, two Array Express datasets, and two private de-identified datasets of patients with systemic inflammatory response syndrome (infection related) from two prospective observational studies—more information about these studies is included in the supplemental material (eSummary 1) (9, 15–24).

Patients were included if

they were adults diagnosed with sepsis or septic shock (regardless of the etiology), and

as gene expression from different RNA extraction kits could introduce bias that may impact confidence in the results, we only included patients whose RNA were extracted using PAXgene Blood RNA tubes (within 24 h of the diagnosis), as it is FDA approved.

Publicly available data used herein is exempt from Institutional Review Board (IRB) review. The private datasets received IRB approval (LSUHSC #8964 and #10074) and written informed consent was obtained from all subjects.

Predictors

Only gene level expression measured across micro-array platforms was used in model development.

Data preparation

Training and validation sets were sourced from thirteen datasets and were constructed using a similar approach to a previous study (25). Datasets which did not included clinical data were used as the training set. Due to the heterogeneous nature of the expression data—and blinded to patient characteristics, treatment, and outcome—custom processing was performed to ensure consistency before modeling. More details are available in the supplemental material (eSummary 2).

Outcomes

The primary outcome was 28-day mortality.

The secondary outcome was the immunological response to hydrocortisone therapy.

Score development and validation

Unsupervised clustering to identify subclasses

We first reproduced previously established subclassification of patients (25). For optimal semantic definition and convenience, the subclasses A, I, and C from the original study were respectively renamed IA (immune-adaptive), IN (immune-innate), and IC (immune-coagulant). The nature and functionality of each subclass was confirmed using Gene Ontology (GO) analysis on both train and validation set (26). Later, the clinical characteristics of the patients were evaluated in the validation dataset to assess their corresponding immune pathophysiology. More information about unsupervised clustering and GO analysis is included in the supplemental material (eSummary 3 and 4).

Generating scores representing the activity level of the immune systems

Once the established subclasses were established via unsupervised clustering, all genes in the training data were extracted and quantile normalization was applied. Then, a set of significantly Differentially Expressed Genes (DEGs) among the previously labelled subclasses (i.e., clusters) was identified. These genes were subsequently down-selected by a forward-search method and optimized by a backward-search method. This process generated minimal sets of up and down predictor genes that could identify patients for each subclass. We then calculated continuous scores for immune-adaptive (IA), immune-innate (IN), and immune-coagulant (IC) using the difference between the geometric mean of the respectively selected up and down predictor genes (27, 28). To confirm the underlying biological meaning of the scores generated, we used the Pearson and Spearman tests to investigate their correlation to the activity level of each component of the immune response, using the complete set of GO genes as the measure of activity. More detail on normalization, DEGs, down-selection, and score correlation analysis is available in the supplemental material (eSummary 5, 6, and 7).

Model performance

To assess if each of the gene-expression scores accurately reflected their respective subclasses, we evaluated the performance of a multi-class classifier (Support Vector Machine) targeting the subclasses by utilizing the scores for overall accuracy and Area Under the Receiver Operator Characteristic Curve (AUROC) using Leave-One-Out cross-validation in the training set. The AUROC was classified as excellent (≥0.9), very good (≥0.8), good (≥0.7), fair (≥0.6), and poor (<0.6).

Immune states

In contrast to the adaptive and innate responses, there is no known effect of hydrocortisone in the coagulation cascade. Consequently, even though the clustering correctly identifies some patients as having a coagulative profile, this classification might not be relevant for the goal of differentiating treatment response to hydrocortisone. Therefore, a better approach could be to identify which of the adaptive versus innate system is prevailing in the patient's physiology. To capture interactions between these two immune systems, we stratified patients into two groups based on their balance using their respective gene-expression scores as a surrogate for their activity levels. Patients showing high adaptive and low innate immune activity were considered immune-adaptive prevalent (IA-P) and patients with low adaptive and high innate immune activity were considered immune-innate prevalent (IN-P). To determine high versus low activity, the median difference between the IA and IN scores was used. Hence, patients with a difference (IA minus IN) higher than the population median were classified as IA-P and otherwise as IN-P. We then assessed the differential of treatment effect in the 28-day mortality of hydrocortisone according to these two groups.

Statistical analysis

For the primary outcome, we assessed the differential treatment response for 28-day mortality according to the interaction between hydrocortisone treatment and the assigned immune-adaptive state (i.e., IA-Prevalent or IN-Prevalent). To do so, we separately analyzed two other independent datasets of patients with vasodilatory shock: Effect of Early Vasopressin vs Norepinephrine as Initial Therapy in Septic Shock (VANISH) trial and Low-dose hydrocortisone reduces norepinephrine duration in severe burn patients: a randomized clinical trial (hereafter referred to as “Burn trial”) (7, 29)—more details on these trials are in the supplementary appendix (eSummary 8 and 9).

For the secondary outcome, we used the gene expression data from the Burn trial to assess the immunological response to hydrocortisone therapy. To do so, we compared the gene-expression scores pre- and 24 h post-hydrocortisone treatment using paired t-tests.

Baseline and outcome data were presented according to the assigned immune-adaptive state (i.e., IA-P or IN-P). Continuous variables were presented as medians with their interquartile ranges and categorical variables as total number and percentage. Proportions were compared using Fisher's exact test and continuous variables were compared using the Wilcoxon rank-sum test. The significance level was set at 0.05 and the software R (v.4.0.2) and HypaHub Big Data Analytics Platform were used for all analyses.

RESULTS

Participants

We included 13 studies from eight different countries in the initial step of our analysis (eTable 1). Of 596 patients that fulfilled our inclusion criteria, 165 were included in the training set, and 431 were included in the validation cohort (eTable 1). Baseline characteristics of the patients in the validation sets are presented in eTable 2. Compared to survivors, non-survivors were older, had more coagulation disorders, higher rate of shock, higher Sequential Organ Failure Assessment (SOFA) scores, and higher Acute Physiology and Chronic Health Evaluation (APACHE) II scores.

Immune-pathophysiological subclasses

A stable optimum of three clusters was the most frequent across the tested space, indicating biological separation at three clusters using 500 genes as the optimal clustering assignment (eFigure 1). One hundred and fifteen samples in the training set were conclusively assigned to one of the groups and the remaining 50 were labeled as unclassified. Gene Ontology (GO) analysis further confirmed that the clusters had distinct biological characteristics on both train and validation sets and were representative of different immune-pathophysiological subclasses (eSummary 10, eFigure 2, eFigure 3).

Analysis of clinical characteristics for the subclasses was done in the validation set using individual-level data. There was no difference between the subclasses in terms of age, gender, and vital signs. However, patients assigned to the immune-innate (IN) subclass had higher white blood cell count, a higher rate of shock, and higher mortality. Patients assigned to the immune-coagulant (IC) subclass had higher APACHE II scores, and a higher rate of coagulation disorders. Immune-adaptive (IA) patients had the lowest white blood cell count, shock rate, APACHE II score, and mortality (eTable 3). A sensitivity analysis of the interactions between APACHE II, shock, and age with immune-adaptive as the reference class confirmed that the subclasses did not simply represent disease severity (immune-innate, P = 0.757; immune-coagulant, P = 0.438).

Gene-expression scores

The detection of DEGs and subsequent application of forward- and backward-search found a set of 15 distinct genes that could identify patients in each subclass. These genes were used to create the IA, IN, and IC scores (eTable 4). Five up-regulated genes (ZNF831, CD3G, MME, BTN3A2, and HLA-DPA1) and one down-regulated gene (STOM) had an AUROC of 0.983 for classifying patients in the IA subclass. For IN, two up-regulated genes (HK3 and SERPINB1) and three down-regulated genes (EPB42, GSPT1, and LAT) had an AUROC of 0.989. Lastly, for IC three up-regulated genes (SLC1A5, IGF2BP2, and ANXA3) and one down-regulated gene (GBP2) had an AUROC of 0.990 (eFigure 4 and eTable5). These results confirmed that gene-expression scores generated from just 15 genes had excellent performance for subclass identification.

Each of the scores were highly representative of the activity level of their assigned immune state according to the respective GO genes: the IA score had a significant correlation with adaptive immune activity (R = 0.77); the IN score had a significant correlation with inflammation activity (R = 0.73); and the IC score had a significant correlation with coagulation activity (R = 0.65) (Fig. 1A–C). The unclassified samples also demonstrated high correlation to one of the immune activity levels (eFigure 5). These additional results demonstrated the scores not only had excellent performance on subclass classification, but also confirm they could be used to represent immune activity levels.

Fig. 1

Scores correlate with gene expression that are related to host-response-state GO terms. (A–C) Figures show that scores by classifier genes have high correlations with score by genes belonging to host-response immune GO terms. The colors represent the subclasses: IA in blue, IN in yellow, IC in grey, and unclassified (U) samples in red.

Immune states

The clinical characteristics of patients classified as immune adaptive prevalent (IA-P) and immune innate prevalent (IN-P) are shown in Table 1. Compared to IN-P patients, IA-P patients had lower white blood cell counts, twice the median lymphocyte count and lower rate of lymphopenia, a lower shock rate, SOFA score, and APACHE II score. Their mortality rate was also lower (22.3% vs. 35.2%, P = 0.004).

Table 1

Validation set clinical features grouped according to immune-adaptive and immune innate prevalent states

Variable	Overall N:431	Immune adaptive prevalent N: 215	Immune innate prevalent N:216	P
Age (years), median [IQR]	61.0 [49.0–74.0]	61.0 [47.0–71.5]	62.0 [50.0–75.5]	0.125
Male, n (%)	240 (55.8)	119 (55.3)	121 (56.3)	0.923
Heart rate (bpm), mean (SD)	109.8 (24.7)	109.4 (23.9)	110.4 (25.8)	0.785
Temperature (oC), median [IQR]	37.7 [36.8–38.8]	37.5 [36.7–38.6]	37.8 [37.1–38.9]	0.177
MAP (mmHg), mean (SD)	75.6 (18.0)	78.3 (17.7)	72.4 (17.9)	0.026
WBC (×103), median [IQR]	14.2 [9.4–20.7]	12.9 [9.0–18.4]	15.3 [9.9–22.9]	0.037
Lymphocytes (cells), median [IQR]	950.0 [516.0–1587.0]	1195.3 [681.8–1781.5]	600.0 [357.2–1292.5]	<0.001
Lymphopenia, n (%)	76 (48.1)	30 (34.9)	46 (63.9)	0.001
Platelets (×103), median [IQR]	189.0 [123.0–258.0]	198.0 [131.0–288.0]	179.0 [102.5–232.5]	0.030
INR, median [IQR]	1.3 [1.2–1.6]	1.3 [1.1–1.5]	1.4 [1.2–1.7]	0.054
Coagulation disorder, n (%)	17 (23.6)	7 (16.7)	10 (33.3)	0.174
Shock, n (%)	197 (54.0)	85 (45.2)	112 (63.3)	0.001
SOFA, mean (SD)	6.2 (3.2)	5.2 (2.8)	7.3 (3.3)	<0.001
APACHE II, mean (SD)	22.1 (8.8)	20.4 (8.9)	24.1 (8.2)	<0.001
Count of Organ Failures, median [IQR]	2.0 [0.5–2.0]	1.0 [0.0–1.0]	2.0 [1.0–3.0]	0.001
Hospital LOS (days), median [IQR]	8.5 [4.2 0–15.0]	8.0 [4.0–12.0]	11.0 [5.0–18.0]	0.021
Non-survivor, n (%)	124 (28.8)	48 (22.3)	76 (35.2)	0.004

APACHE II score, Acute Physiology and Chronic Health Evaluation; coagulation disorder, Platelet < 100,000 and INR > 1.3; IQR, Interquartile range; LOS, length of stay; lymphopenia, lymphocytes ≤ 900 cells; MAP, Mean arterial blood pressure; SD, standard deviation; SOFA score, Sequential Organ Failure Assessment; WBC, white blood cell.

Differential treatment response

Hydrocortisone was associated with higher mortality rate only for IA-P patients in both studies: significantly in the VANISH (43.3% vs. 14.7, P = 0.028) and, numerically in the Burn trial (57.1% vs. 0.0%, P = 0.99) as shown in Table 2.

Table 2

Differential treatment response on VANISH and Burn trials

	Immune-adaptive prevalent (IA-P)		Immune-innate prevalent (IN-P)
	Hydrocortisone	Placebo	Hydrocortisone	Placebo	Interaction P
VANISH TRIAL
n	30	34	28	25
28-day mortality (%)	43.3	14.7	32.1	40	0.028
BURN trial
n	7	8	8	7
28-day mortality (%)	57.1	0.0	25	28.6	0.99

Immunological response to hydrocortisone

To assess the immunobiological response to hydrocortisone, we analyzed the pre- and 24 h post-treatments expression profiles of the patients included in the Burn trial. Twenty-four hours post-intervention, hydrocortisone led to a significant decrease of the immune adaptive score (P = 0.033) as shown in Figure 2, leading to a drop in the percentage of immune adaptive prevalent patients compared to placebo (eFigure 6).

Fig. 2

Classifier scores can be used to classify patients as immune-adaptive prevalent (IA-P) and immune-innate prevalent (IN-P) in the Burn trial. Patients were stratified to IA-P and IN-P groups based on baseline samples. The line marked by an asterisk shows that hydrocortisone treatment leads to a significant decrease of immune adaptive score of IA-P patients comparing baseline samples and 24 h post-treatment samples but not in placebo groups (P = 0.033). The blue and yellow lines represent the hydrocortisone and placebo groups, respectively.

These observations were further confirmed by DEG analysis between pre- and 24 h post-treatment. DEG results showed that the hydrocortisone group had more gene expression changes compared to the placebo group (Fig. 3). Additionally, GO analysis revealed that the genes down-regulated by hydrocortisone related to adaptive immunity pathways, such as T-cell activation and antigen receptor-mediated signaling, an effect that was not observed in the placebo group (Fig. 3).

Fig. 3

Gene Ontology (GO) analysis on differentially expressed genes between pre- and 24 h post-treatment. Enrichment of Gene Ontology terms in significantly up-regulated and down-regulated genes (P-value < 0.01 and fold change > 1.1) 24 h after the intervention (i.e., Hydrocortisone or placebo). The size of dot represents the enrichment odds ratio of GO terms. The color represents the adjusted P-value of GO term enrichment. The plot indicates hydrocortisone group has more significant DEGs between pre- and 24 h post-treatment than placebo group. Additionally, the GO enrichment analyses infer that hydrocortisone treatment might uniquely lead to expression decrease of genes involved in adaptive immunity pathways.

DISCUSSION

Key findings

We applied unsupervised clustering of gene expression data to identify biological subclasses of vasodilatory shock. We found that continuous gene-expression scores generated using 15 genes can accurately reflect the activity level of immune states. In patients with vasodilatory shock caused by either sepsis or severe burns, these scores identified those with a prevalent immune-adaptive state. Such patients from the VANISH trial had a significantly greater 28-day mortality when treated with hydrocortisone instead of placebo and similar effect in the Burn trial. Moreover, compared to baseline, after 24 h, the immune-adaptive (IA) score of patients treated with hydrocortisone decreased, but remained unchanged with placebo.

Relationship to previous findings

Whole blood RNA expression profiling of sepsis patients is widely used to discover sepsis subtypes since whole blood consists of cells associated with the innate and adaptive immune systems which play critical roles in sepsis. Three other groups have found the same differential treatment effect in the VANISH trial (9, 30, 31). It is compelling to note that all these studies used different datasets, populations and clustering techniques, corroborating the robustness of our findings. However, none of these previous studies demonstrated their model could accurately reflect immune activity levels, nor provide the impact of the hydrocortisone on disease progression using the patient's gene expression. In addition, Sinha et al. have identified two COVID-19-related ARDS subgroups with differential treatment response to corticosteroids which further supports our findings (32).

The use of continuous scores to find responders to therapies is not new in the medical field. In oncology, investigators have already demonstrated that immune related antigen scores are able to identify patients with better response rates to immune therapies (33, 34). Experimental factors such as sample preparation and human manual errors can introduce variability to the gene expression. To deal with this problem, our gene-expression score was created using the geometric mean of up and down predictor genes. Both the use of geometric mean and the differential expression technique have been widely used to reduce noise (27, 28), potentially leading to a better representation of the trends in gene expression for each subclass.

Implications of study findings

It is widely accepted that the immune system is tightly regulated with both adaptive and innate states as major components. Our findings imply that assessment of such immune states can now be achieved with gene-expression scores. They also suggest that hydrocortisone may harm patients with a predominant immune-adaptive state, an effect consistent with the adaptive immunosuppressive properties of hydrocortisone (35). Finally, our findings suggest that gene expression scores may provide the means to measure the activity level of immune system components and new biological insights toward the development of personalized treatment.

Study strengths and limitations

This study has several strengths. To the best of our knowledge, this is the first study using continuous gene expression scores to show differential treatment response in patients with vasodilatory shock. To do so, we used data from two rigorous randomized controlled clinical trials, both assessing hydrocortisone therapy versus placebo with paired transcriptomic data collected close to the time of randomization. Additionally, we also show how the immune adaptive score, which is representative of the adaptive immune system's state, is impacted by the natural progression of the inflammatory response syndrome in vasodilatory shock patients, both with and without hydrocortisone intervention. This provides a potential biological explanation of our findings. Moreover, this study introduces a new concept—the gene-expression scores—that may help the research community to evaluate new and previously attempted therapeutic approaches in more homogenous groups of patients with respect to their immunological characteristics. Lastly, we foresee that the use of a real time reverse transcription polymerase chain reaction (RT-qPCR) assay will allow a peripheral whole blood specimen to be collected from patients and added directly to a cartridge that performs RNA extraction. This process will allow measurement of the gene expression for the 15 key genes described here in approximately 1 h and enable the application of personalized medicine in the care of vasodilatory shock patients.

We acknowledge several limitations. First, we have included only the samples that lack paired clinical data in the training set. However, this shortcoming is unlikely to have impacted the results of our study. This is because our unsupervised model that generated the sepsis subclasses included only gene expression data and used a large and complete clinical validation set. Second, due to a limitation in the COCONUT framework that requires healthy samples to be present in the dataset to perform normalization, we could not directly assess the unsupervised clustering in the validation set. Fortunately, this limitation is mitigated using leave-one-out sampling when generating the AUROC curves, as well as by the GO and clinical analyses performed in the validation set. Third, the Burn Trial was interrupted before the original sample size was obtained, which limited our analysis. Fourth, the number of patients assessed in the VANISH and especially Burn trials is small and may have exposed our findings to the risk of type I error. Fifth, we assumed that all vasodilatory shock patients included in this analysis had a high cardiac output as it is a dominant physiologic characteristic of this syndrome. However, this assumption cannot be verified and is a limitation of our analysis. Sixth, the VANISH trial data is further limited by a subset of patients who were on high doses of vasopressors and thus eligible for hydrocortisone therapy. Therefore, this study is hypothesis generating due to the limited generalizability of our findings to a broader population of patients with vasodilatory shock. Finally, these findings were derived from a secondary analysis of past trials and thus prospective validation is still required.

CONCLUSIONS

The findings of this study suggest that a gene-expression score derived from sepsis subclasses may be used to guide therapy in patients with vasodilatory shock. In particular, patients with a prevalent immune-adaptive state may be harmed by treatment with hydrocortisone. This observation implies that the use of gene-expression scores may enable researchers to measure the activity level of immune system components and provide new biological insights toward the development of personalized treatment.