A novel immuno-oncology algorithm measuring tumor microenvironment to predict response to immunotherapies.

Immune checkpoint inhibitor (ICI) therapies can improve clinical outcomes for patients with solid tumors, but relatively few patients respond. Because ICI therapies support an adaptive immune response, patients with an active tumor microenvironment (TME) may be more likely to respond, and thus biomarkers capable of discerning an active from a quiescent TME may be useful in patient selection. We developed an algorithm optimized for genes expressed in the mesenchymal and immunomodulatory subtypes of a 101-gene triple negative breast cancer model (Ring, BMC Cancer, 2016, 16:143) as a means to capture the immunological state of the TME. We compared the outcome of the algorithm (IO score) with the 101-gene model and found 88% concordance, indicating the models are correlated but not identical, and may be measuring different TME features. We found 92.5% correlation between IO scores of matched tumor epithelial and adjacent stromal tissues, indicating the IO score is not specific to these tissues, but reflects the TME as a whole. We observed a significant difference in IO score (p = 0.0092) between samples with high tumor-infiltrating lymphocytes and samples with increased neutrophil load, demonstrating agreement between IO score and these two prognostic markers. Finally, among non-small cell lung cancer patients receiving immunotherapy, we observed a significant difference in IO score (p = 0.0035) between responders and non-responders, and a significant odds ratio (OR = 5.76, 95% CI 1.30-25.51, p = 0.021), indicating the IO score can predict patient response. The immuno-oncology algorithm may offer independent and incremental predictive value over current biomarkers in the clinic.

Introduction

Immunotherapies using immune checkpoint inhibitors (ICIs) are now standard of care in the treatment of lung cancer, breast cancer, and other solid tumor types [1, 2]. Although ICIs are able to improve clinical outcomes for patients with a variety of solid tumors, only a small subset of patients respond [3, 4]. Moreover, ICIs can cause immune-related adverse events, some of which are clinically serious and potentially life-threatening [5, 6]. Hence, the ability to identify patients who are more likely to benefit from ICI therapy with minimal toxicity is essential for optimizing immunotherapy use.

Predictive biomarkers that can support patient selection for ICI therapy are limited. Expression of programmed death-ligand 1 (PD-L1) on tumor cells has been most widely investigated, but PD-L1 testing does not consistently predict patient benefit from immunotherapy [7, 8, 9]. The lack of effective biomarkers for predicting response to ICIs is partially due to the complexity of the tumor-immune system interactions of the tumor microenvironment (TME). Within the TME is a complex and dynamic milieu of non-malignant cells that interact with each other and with the tumor cells, affecting tumor growth, invasion and metastasis [10, 11]. As such, biomarkers that are able to capture the complex interactions of the TME could be more useful in selecting patients who will benefit from ICI therapies.

In previous work, a 101-gene model was developed to clinically subtype triple negative breast cancer (TNBC) patients [12]. This algorithm classified TNBC into five molecular subtypes, including two basal like (BL1 and BL2), luminal androgen receptor (LAR), mesenchymal (M), and mesenchymal stem-like (MSL), with each of these subtypes further classified by a positive or negative immunomodulatory (IM) signature. In applying the 101-gene model to independent TNBC cohorts, we noticed that tumors of the M subtype never had a positive IM signature, an observation that was consistent with studies showing that the M and IM subtypes are inversely correlated [13, 14]. This observation led to the hypothesis that the M subtype would be antithetical to the IM subtype, with the former indicating a more quiescent immunological state and the latter indicating an immunologically active state. Additionally, we hypothesized that the molecular basis for the M and IM subtypes would translate across other solid tumor types based on features of the TME that are driving this profile. We reasoned that, if this were true, it may be possible to develop a gene expression algorithm to measure the TME by optimizing a gene set to include those most relevant to the M and IM subtypes. Such an algorithm could potentially distinguish tumors in a more quiescent state from those that are immunologically active and may have expanded utility across multiple cancer types.

Here we describe the development and validation of a novel immuno-oncology algorithm that measures the TME. This algorithm was optimized using genes expressed in both quiescent and immunologically active tumors and may be useful in predicting response to immunotherapies.

Material and methods

Data analysis

All analyses, unless otherwise stated, were done on RStudio Version 1.2 utilizing R version 3.6 [15,16].

Gene expression dataset processing

Twenty-five gene expression profile data sets (Table 1), representing three microarray platforms, were downloaded from the publicly available Gene Expression Omnibus (GEO, ncbi.nlm.nih.gov/geo/). Data were combined from raw microarray expression (CEL) files collectively normalized by robust multiarray average (RMA), and log transformed. Samples from this data set were pared down to triple negative status using a bimodal distribution of ESR1, ERBB2, and PGR genes, resulting in 1284 unique TNBC samples. Of these, 994 unique TNBC samples were used to train the model, and the remaining 335 unique TNBC samples were used for model validation.

Table 1

Source of TNBC specimens for model Training and Validation.

Dataset	TNBC Specimens	Dataset	TNBC Specimens
GSE1456	44	GSE7904	17
GSE1561	21	GSE10780	5
GSE2034	59	GSE11121	21
GSE2109	55	GSE12093	57
GSE2603	35	GSE12763	5
GSE2990	11	GSE13787	10
GSE3494	27	GSE16716	62
GSE3744	17	GSE25066	178
GSE5327	35	GSE31519	67
GSE5364	36	GSE58812	107
GSE5462	2	GSE76124	198
GSE6596	8	GSE76250	165
GSE7390	42

For genes represented by multiple probes, the probe with the highest inter-quartile range was selected to prioritize genes with a large dynamic range of expression. Batch correction was performed using an Empirical Bayes method, ComBat [17]. Patient datasets were previously made publicly available under the ethical policies of the National Institutes of Health's Gene Expression Omnibus (GEO) database. No additional ethics review was required for the in-silico analysis of these datasets.

Model building

Model building for the novel immuno-oncology algorithm was performed using R version 3.5.2 (Figure 1). As suggested by Ring et al., large scale gene expression classifiers have an inherent problem with the inclusion of genes which have minimal signal contribution across each classifier [12]. Therefore, a further reduction of the genes from the 101-gene signature to enrich for those which most strongly contribute to the immunomodulatory and mesenchymal classifiers would hypothetically lead to a more clinically robust and practical biomarker.

Figure 1

Gene selection process for building the novel immuno-oncology algorithm. Gene set resulted from data set normalization, batch correction, gene set enrichment analysis, and elastic net modeling.

To this end, the 101-gene signature was used to identify gene sets that distinguished the classes via gene set enrichment analysis (GSEA) using the C2 curated gene sets of canonical pathways [18]. Elastic-net regularized linear models were employed to create individual subclassifying models for the BL1, BL2, LAR, MSL, M, and IM subtypes with the subtypes treated as a multinomial variable [19]. The 30 genes utilized for the M and IM subtype classifications with this model were then used to derive a logistic elastic net model on the new data set, minus three genes whose probes had been reassigned between analyses (UBD, JCHAIN, and KRT17). Strength of association with classification variables was assessed using ten-fold cross validation of the misclassification error. The model threshold for determining the immuno-oncology score (IO score) was determined using the area under the curve (AUC) [20] to maximize sensitivity to the immunomodulatory score from 101-gene signature, in contrast to the significance of the correlation method for determining threshold previously described by Ring et al. [12].

GSE81838 dataset analysis of TNBC tumor epithelial and adjacent stromal tissue

Microarray data were obtained from GSE81838 where laser-capture microdissection had been performed on 10 TNBC tumors to isolate malignant epithelial cell-enriched areas and the adjacent stromal cell-containing areas of the tumor sections [13]. The IO scores for each sample were obtained and correlated between the matched tumor epithelial and adjacent stromal tissue using Spearman's method.

TCGA breast cancer datasets and analysis

Gene expression profiles from breast cancer specimens collected for The Cancer Genome Atlas (TCGA) were obtained from the National Cancer Institute Genomic Data Commons Data Portal. TNBC status was confirmed by bimodal modeling of ESR1, PGR, and ERBB2 gene expression, resulting in 180 total samples with matching tumor infiltrating lymphocytes (TILs) presence and intensity as described previously [13]. Neutrophil presence was obtained by the TCGA study investigators who assessed the hematoxylin and eosin (H&E) staining to discern neutrophil presence and was aligned to the TNBC samples. The IO scores of samples with intense TIL staining and samples with neutrophil presence of 30% or greater was assessed by the Welch t-test for significance.

GEO non-small cell lung cancer (NSCLC) datasets and analysis

The clinical response to anti-PD-1 therapy and expression data of advanced NSCLC patients in the GSE135222 (27 patients) and GSE126044 (16 patients) cohorts were obtained from GEO. Response was measured in both cohorts using Response Evaluation Criteria in Solid Tumors (RESCIST) metrics [21], where patients exhibiting partial response or stable disease for >6 months were classified as responders [22, 23]. Because response was defined in the same manner for both cohorts, we were able to combine the data for purposes of the analysis. Expression data from the combined cohort were processed using the novel immuno-oncology algorithm and analyzed by IO score. The difference in IO score between responders and non-responders was evaluated for significance using the Welch t-test. The data from the combined cohort was then evaluated for the correlation of IO score to objective response. The predefined threshold was used to divide patients into IO score positive and negative and compared to objective response to calculate an odds ratio.

Results

Concordance between IO score from novel algorithm and IM status from 101-gene model

We hypothesized using features of the M and IM subtypes from the 101-gene model could be combined in an independent algorithm with a reduced gene set that may increase the likelihood of capturing biological processes relevant to immunotherapy response. To this end, an independent expression-based centroid model, defined by the M and IM features of the 101-gene model, was obtained from elastic net modeling and contained a total of 27-genes (Table 2).

Table 2

IO signature 27-gene list.

HGNC_Gene_symbol	ensmbl_ID
APOD	ENSG00000189058
ASPN	ENSG00000106819
CCL5	ENSG00000271503
CD52	ENSG00000169442
COL2A1	ENSG00000139219
CXCL11	ENSG00000169248
CXCL13	ENSG00000156234
DUSP5	ENSG00000138166
FOXC1	ENSG00000054598
GZMB	ENSG00000100453
HTRA1	ENSG00000166033
IDO1	ENSG00000131203
IL23A	ENSG00000110944
ITM2A	ENSG00000078596
KMO	ENSG00000117009
KRT16	ENSG00000186832
KYNU	ENSG00000115919
MIA	ENSG00000261857
PSMB9	ENSG00000240065
PTGDS	ENSG00000107317
PLAAT4	ENSG00000133321
RTP4	ENSG00000136514
S100A8	ENSG00000143546
SFRP1	ENSG00000104332
SPTLC2	ENSG00000100596
TNFAIP8	ENSG00000145779
TNFSF10	ENSG00000121858

We then compared the independent immuno-oncology algorithm and 101-gene model by testing a validation cohort of 335 unique TNBC samples (Table 3). The resulting matrix indicates 88% concordance for IO+/IM+ and IO-/IM- scores, indicating that the novel immuno-oncology algorithm and 101-gene model are strongly correlated but not identical, and may be measuring different features of the TME.

Table 3

Concordance between IM status from the 101-gene model and IO score from the novel immuno-oncology algorithm within the validation cohort of 335 unique TNBC samples.

		101-gene TNBC Model
		IM+	IM-
IO Algorithm	IO+	82 (24%)	37 (11%)
	IO-	2 (1%)	214 (64%)

Correlation of IO score to tumor epithelial and adjacent stromal tissue in TNBC

The IO score describes the TME which generally includes both tumor and stromal tissue. As such, we hypothesized that both tumor and stromal tissue would contribute to the IO score in describing the state of the TME. As a test of this hypothesis, we analyzed the IO scores from matched TNBC tumor epithelial and adjacent stromal tissue samples using the GSE81838 dataset. This dataset includes gene expression data from 10 patients with TNBC tumors that had been subjected to laser-capture microdissection to parse the epithelial tumor mass from the adjacent stromal tissue resulting in 20 unique samples. Due to the limited sample size (only 20 samples) and replicate patient identification (only 10 patients), we chose to evaluate concordance of the IO score of matched tumor epithelial and adjacent stromal tissue samples for each patient using Spearman's method. We observed a 92.7% (p < 0.001) correlation of IO scores between the tissue types when matched to each unique patient, suggesting that the IO score is a measure of the tumor microenvironment spanning the tumor and stromal regions, and may tolerate variable tumor fraction in available samples to obtain a reliable signal.

Correlation of IO score to high levels of TILs or neutrophils in TNBC

The presence of high levels of TILs indicates an active immunological state and is considered prognostic of good outcome with ICIs in TNBC patients [24]. In contrast, increased neutrophil load in the tumor may be contraindicative of response to immunotherapy [25, 26]. Therefore, as a biological test of whether the IO score may indicate an active immunological state or not, we evaluated the IO scores of TNBC samples with high TILs and samples with increased neutrophil load obtained from TCGA (Figure 2). Using the Welch t-test, we found that the difference in IO score between TNBC samples with high TILs (0.09) and samples with increased neutrophil load (-0.30) was significant (p = 0.0092), indicating agreement with these prognostic markers. These results suggest that a positive IO score may possess features associated with a favorable outcome to immunotherapy, while a negative IO score may indicate poor immunotherapy response.

Figure 2

Box and Whisker plot displaying IO scores of TNBC samples from TCGA with high levels of TILs as compared to samples with increased neutrophil load. The IO score threshold is indicated at 0.09. The line within the box plots represents the median and the cross represents the mean.

Correlation of IO score to immunotherapy response in NSCLC patients

As a test of whether the IO score may hold prognostic potential for other sold tumor types, we evaluated the IO scores from a combined cohort of NSCLC patients where response was defined as exhibiting partial response or stable disease for >6 months. The average IO score of responders (0.29) and non-responders (-0.096) was found to be significantly different by the Welch t-test (p = 0.0035) (Figure 3), clearly separated on both sides of the predefined threshold. Notably, while the immuno-oncology algorithm was generated using RNA expression data from microarrays, the input data for these NSCLC cohorts were generated using whole transcriptome RNA-Seq, demonstrating the ability of the immuno-oncology algorithm to analyze data from multiple platforms.

Figure 3

Box and Whisker plot displaying IO scores from Responders (R) and Non-Responders (Non-R) from the combined NSCLC cohorts. The IO score threshold is indicated at 0.09. The line within the box plots represents the median and the cross represents the mean.

We further evaluated the prognostic potential of the IO score in NSCLC by calculating the Odds Ratio (OR), finding that clinical response to immunotherapy was significantly associated with a positive IO score (OR = 5.76, 95% CI 1.30 to 25.51, p = 0.021) using the predefined threshold for IO positivity.

Discussion

Immunotherapy with ICIs has revolutionized cancer treatment and is responsible for a significant portion of the largest single-year drop in cancer mortality rates from 2016 to 2017 [27]. Despite its effectiveness in multiple cancer types, only a small proportion of patients respond to ICI therapy. Recent estimates of U.S. cancer patients eligible for ICI therapy have been lowered from 12.5% to 10.9%, in part due to the failure of ICI drugs to show improvement in overall survival or progression-free survival [28]. In order to realize the full potential of ICI therapies, reliable biomarkers are necessary to distinguish patients who will benefit and to avoid unintended morbidities.

Capturing the immunological state of the TME appears to be a critical pillar of predicting response to immunotherapies [29]. Whereas the active TME is defined as having an established innate and adaptive immune response to the tumor as well as supportive immune recruitment by the surrounding stroma, the quiescent state contains more immunosuppressive features including increased extracellular matrix components [30] and an increased pro-inflammatory response enhancing vascularization [31] which plays a tumor supportive role. Because ICI therapies typically support an adaptive immune response, patients with an active TME may be more likely to respond [32]. Current gold standard biomarkers only capture a defined aspect of the TME and have limited ability to predict immunotherapy response. Thus, there is an unmet clinical need for biomarkers capable of discerning an active from a quiescent TME which may be useful in predicting response to ICI therapy.

Here we have described a novel immuno-oncology algorithm that was optimized for genes expressed in the M and IM subtypes of TNBC as a means to capture the immunological state of the TME (Figure 4). Earlier studies have shown that the IM subtype may have predictive value in immunotherapies, and that the M and IM subtypes are inversely correlated [13, 14]. We hypothesized that the IM subtype reflects an immunologically active TME which is indicative of response to immunotherapies, whereas the M subtype reflects a more quiescent TME. We reasoned that the M and IM subtypes could inform an algorithm to measure the TME immunological state that may be useful in selecting patients likely to benefit from ICI therapies.

Figure 4

Overview of IO score as a measure of the quiescent or immunologically active state of the tumor microenvironment (TME). We hypothesized that a negative IO score may indicate a quiescent state, where the tumor cells are more actively promoting angiogenesis, inducing an inflammatory response, and stimulating cancer-associated fibroblasts which collectively is constructing extracellular matrix. By comparison, a positive IO score may indicate an immunologically active TME with reduced inflammatory characteristics combined with an increase in the innate and adaptive immune systems increasing tumor cell invasion. Whereas a biomarker such as PD-L1 may be present in both states, the immuno-oncology algorithm is able to distinguish a quiescent from an active TME.

In developing the immuno-oncology algorithm, we adapted some of the methods used to derive the 101-gene model described by Ring et al. to generate shrunken centroids [12]. As hypothesized by Ring, feature and dimensional reduction of large-scale gene expression datasets minimizes overfitting which can improve clinical subtype identification. This hypothesis seems to be supported by the outcome of the validation cohort in which we observed 88% concordance between the IM score of the 101-gene model and the IO score of the immuno-oncology algorithm. Among the discordant group, 11% is represented by samples that were IO positive but IM negative. Further studies are planned to evaluate the performance of the novel immuno-oncology algorithm in clinical cohorts.

Several observations support our hypothesis that the novel immuno-oncology algorithm measures the immunological state of the TME. For one, using microarray data from GSE81838 we were able to obtain a 92.5% Spearman correlation coefficient between IO scores of the tumor epithelial and adjacent stromal tissues. These results indicate that the IO score is not specific to either of these tissue types, but rather reflects the broader TME as a whole. In another line of experiments with TCGA breast cancer datasets, we observed a significant difference in IO score (p = 0.0092) between TNBC samples with high TILs and samples with increased neutrophil load. This result demonstrates agreement between IO score and these two prognostic markers, and suggests that a negative IO score possesses features consistent with poor response to immunotherapy, whereas a positive IO score correlates with features associated with a favorable outcome.

To further elucidate how the different 27-genes capture a complex signature of the TME and simplify to an indicator of response to immunotherapies we have separated the genes by whether they impact immunomodulatory (hot) or mesenchymal (cold) signatures of the IO score (Table 4). While CDH1 (E-cadherin) and TGFB1 (TGF-β) are commonly associated with the epithelial to mesenchymal transition (EMT), they are not specifically included in our algorithm. However, many of the cold genes of the 27-gene algorithm appear to interact with the EMT pathway suggesting this is an important mechanism we are capturing.

Table 4

Description of 27-gene IO algorithm gene profile.

‘Hot’ TME genes	CCL5, CD52, CXCL11, CXCL13, GZMB, IDO1, IL23A, ITM2A, KMO, KYNU, PSMB9, PTGDS, RARRES3, RTP4, S100A8, TNFAIP8, TNFSF10
‘Cold’ TME genes	APOD, ASPN, COL2A1, FOXC1, HTRA1, KRT16, MIA, SFRP1
Uncertain	DUSP5, SPTLC2

Finally, we were able to demonstrate that the IO score is informative not only in TNBC but also in NSCLC. Using RNA expression data matched to clinical outcome in a combined cohort of NSCLC patients treated with an immunotherapy drug, we observed a significant difference in IO score (p = 0.0035) between responders and non-responders. In addition, we observed a significant odds ratio (OR = 5.76, 95% CI 1.30 to 25.51, p = 0.021), indicating the IO score is capable of predicting response to immunotherapies. Together, these data demonstrate the potential clinical utility of the immuno-oncology algorithm in providing a means to optimize immunotherapy use by predicting response to ICIs. We plan to further evaluate the performance of the immuno-oncology algorithm in clinical cohorts.

Conclusions

There is an unmet need for improved biomarkers to optimize ICI immunotherapy use in clinical settings. The novel immuno-oncology algorithm could potentially address this need by providing a means to distinguish patients likely to benefit from treatment with ICIs. Unlike previously described biomarker models, the novel immuno-oncology algorithm measures the immunological state of the TME as a means to capture the interplay of the patient's immune system and tumor immune evasion. The concept that “tumors are wounds that do not heal” has been used to describe this interplay as the tumor co-ops the wound healing response which encompasses immunosurveillance as well as various aspects of wound healing that appear to be components of tumor maintenance and growth [33]. What makes the immuno-oncology algorithm unique is its ability to capture aspects of immunosurveillance, immunosuppression, and immune evasion as a tumor transitions from a proliferative to a metastatic state. Moreover, the immuno-oncology algorithm is capable of analyzing data obtained from multiple platforms. We believe that by measuring the immunological state of the TME as a whole, the novel immuno-oncology algorithm may offer independent and incremental predictive value over the current gold standard biomarkers in the clinic.

Declarations

Author contribution statement

Tyler J. Nielsen, Brian Z. Ring, Robert S. Seitz: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

David R. Hout: Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Brock L. Schweitzer: Analyzed and interpreted the data; Wrote the paper.

Funding statement

This work was supported by Insight Genetics, Inc. a subsidiary of Oncocyte Corporation.

Data availability statement

Data included in article/supplementary material/referenced in article.

Declaration of interests statement

T.J.N., B.L.S., and D.R.H. are employed by Oncocyte Corporation, the commercial entity that markets the immuno-oncology algorithm as DetermaIO®. B.Z.R. and R.S.S. are consultants contracted by Oncocyte Corporation. We have filed a provisional patent around some of the work cited in this manuscript.

Additional information

No additional information is available for this paper.