key: cord-1037846-q5fr9m76
authors: Yao, Lijing; Rey, Diego Ariel; Bulgarelli, Lucas; Kast, Rachel; Osborn, Jeff; Van Ark, Emily; Fang, Li Tai; Lau, Bayo; Lam, Hugo; Teixeira, Leonardo Maestri; Neto, Ary Serpa; Bellomo, Rinaldo; Deliberato, Rodrigo Octávio
title: Gene Expression Scoring of Immune Activity Levels for Precision Use of Hydrocortisone in Vasodilatory Shock
date: 2022-01-25
journal: Shock
DOI: 10.1097/shk.0000000000001910
sha: b940f9e57147813a5f359929190251df03bbb49a
doc_id: 1037846
cord_uid: q5fr9m76

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.

Vasodilatory shock is the most common form of shock in Intensive Care Unit (ICU) patients (1) (2) (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) (11) (12) (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.

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, http://links.lww.com/SHK/B407) (9, (15) (16) (17) (18) (19) (20) (21) (22) (23) (24) .

Patients were included if 1. they were adults diagnosed with sepsis or septic shock (regardless of the etiology), and 2. 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.

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

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, http://links.lww.com/SHK/B407).

The primary outcome was 28-day mortality. The secondary outcome was the immunological response to hydrocortisone therapy.

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, http://links.lww.com/SHK/B407).

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), immuneinnate (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, http://links.lww.com/SHK/B407).

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 immuneadaptive 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.

For the primary outcome, we assessed the differential treatment response for 28day 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, http://links.lww.com/SHK/B407).

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.

We included 13 studies from eight different countries in the initial step of our analysis (eTable 1, http://links.lww.com/SHK/ B407). 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, http://links.lww.com/SHK/B407). Baseline characteristics of the patients in the validation sets are presented in eTable 2, http://links.lww.com/SHK/B407. 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.

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, http://links.lww.com/SHK/B407). 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, http://links.lww.com/SHK/B407).

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, http://links.lww.com/SHK/B407). 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 (immuneinnate, P ¼ 0.757; immune-coagulant, P ¼ 0.438).

The detection of DEGs and subsequent application of forwardand 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, http://links.lww.com/SHK/ B407). 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 upregulated genes (SLC1A5, IGF2BP2, and ANXA3) and one downregulated gene (GBP2) had an AUROC of 0.990 (eFigure 4 and eTable5, http://links.lww.com/SHK/B407). 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, http://links.lww.com/SHK/B407). 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. 

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).

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 .

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, http://links.lww.com/SHK/B407).

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) .

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 immuneadaptive 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.

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.

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.

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 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. 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.

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.

Comparison of dopamine and norepinephrine in the treatment of shock

Definitions and pathophysiology of vasoplegic shock

Transcriptome modulation by hydrocortisone in severe burn shock: ancillary analysis of a prospective randomized trial

Immune therapy in sepsis: are we ready to try again?

Severe sepsis and septic shock

Low-dose hydrocortisone reduces norepinephrine duration in severe burn patients: a randomized clinical trial

Sepsis: multiple abnormalities, heterogeneous responses, and evolving understanding

Transcriptomic signatures in sepsis and a differential response to steroids. From the VANISH randomized trial

Classification of patients with sepsis according to blood genomic endotype: a prospective cohort study

Genomic landscape of the individual host response and outcomes in sepsis: a prospective cohort study

Robust classification of bacterial and viral infections via integrated host gene expression diagnostics

Developing a clinically feasible personalized medicine approach to pediatric septic shock

Sepsis subclasses: a framework for development and interpretationÃ

Transcriptomic correlates of organ failure extent in sepsis

Aberrant cell cycle and apoptotic changes characterise severe influenza A infection-a meta-analysis of genomic signatures in circulating leukocytes

A distinct influenza infection signature in the blood transcriptome of patients with severe community-acquired pneumonia

Early and dynamic changes in gene expression in septic shock patients: a genome-wide approach

An integrated transcriptome and expressed variant analysis of sepsis survival and death

Genomic transcriptional profiling identifies a candidate blood biomarker signature for the diagnosis of septicemic melioidosis

Prevention of nosocomial infections in critically ill patients with lactoferrin: a randomized, double-blind, placebocontrolled study

Torcetrapib induces aldosterone and cortisol production by an intracellular calcium-mediated mechanism independently of cholesteryl ester transfer protein inhibition

Identifying key regulatory genes in the whole blood of septic patients to monitor underlying immune dysfunctions

Excessive innate immune response and mutant D222G/N in severe A (H1N1) pandemic influenza

Unsupervised analysis of transcriptomics in bacterial sepsis across multiple datasets reveals three robust clusters

ToppGene Suite for gene list enrichment analysis and candidate gene prioritization

Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes

The logarithmic transformation and the geometric mean in reporting experimental IgE results: what are they and when and why to use them?

Effect of early vasopressin vs norepinephrine on kidney failure in patients with septic shock: the VANISH randomized clinical trial

Deep learning-based clustering robustly identified two classes of sepsis with both prognostic and predictive values

External corroboration that corticosteroids may be harmful to septic shock Endotype A patients

Latent class analysis reveals COVID-19-related acute respiratory distress syndrome subgroups with differential responses to corticosteroids

Analysis of 100,000 human cancer genomes reveals the landscape of tumor mutational burden

Analytic validation of tumor mutational burden as a companion diagnostic for combination immunotherapy in non-small cell lung cancer

Corticosteroids are associated with repression of adaptive immunity gene programs in pediatric septic shock

This article is dedicated to the memory of Hector R. Wong MD, who was a great mentor to many investigators including authors of this paper.