A robust panel based on tumour microenvironment genes for prognostic prediction and tailoring therapies in stage I-III colon cancer.

BACKGROUND: Tumour microenvironment (TME) is critical for the regulation of cancer development as well as therapy. The objective of the current study was the development of a robust prognostic model based on TME-relevant genes.
METHODS: Five public microarray datasets providing clinical information were obtained. The least absolute shrinkage and selection operator regression method was used to reduce the dimensionality of robust prognostic genes identified via the bootstrap method.
FINDINGS: We established a prognostic panel, designated as tumour microenvironment risk score (TMRS), consisting of 100 genes. With specific risk score formulae, the TMRS panel possesses a strong ability to predict relapse-free survival and overall survival through both univariate and multivariate analyses. Compared with the TNM stage, the TMRS panel showed much higher predictive accuracy. Further analysis revealed that patients with higher TMRS scores exhibited no therapeutic benefits from adjuvant chemotherapy, probably due to the activation of stromal relevant pathways and infiltration of stromal cells. Besides colon cancer, the TMRS panel was also revealed to be a reliable tool for prognostic prediction and chemotherapeutic decision-making in gastric cancer. Its value in predicting immunotherapy outcomes was also confirmed in two other cohorts consisting of metastatic urothelial carcinoma patients and melanoma patients.
INTERPRETATION: Our TMRS panel may be an effective tool for survival prediction and treatment guidance in patients with stage I-III colon cancer. FUND: This work was supported by the National Natural Science Foundation of China (No. 81772580) and Guangzhou Planed Project of Science and Technology (No. 201803010070).

Evidence before this study

The heterogeneity of tumour microenvironment (TME) contains multiple dimensions of information on patient prognosis and treatment response. Currently, several signatures based on TME genes have been reported. However, as the genes included in these signatures were not all correlated with prognosis, the clinical practicality of these signatures was not satisfactory. In colon cancer, recent advances in high-throughput gene testing technology have led to the development of some molecular signatures for prognosis prediction and personalisation of treatment paradigms. However, to the best of our knowledge, none of these signatures were established based on TME-relevant genes. In addition, an immunoscore system has been developed based on IHC staining of CD3 and CD8. Unfortunately, the accuracy of prognosis prediction using the immunoscore system was found to be limited, as reported by the latest multi-central clinical research. Therefore, it is necessary to.

develop a novel TME-related prognostic model to improve predictive accuracy and identify patients for whom chemo- and immunotherapies may be more beneficial.

Added value of this study

We developed a robust prognostic panel utilising the machine learning method based on TME-relevant genes for stage I–III colon cancer patients, designated as the “tumour microenvironment risk score (TMRS)”. This panel not only accurately predicted relapse-free survival and overall survival among colon cancer patients, but also served as a biomarker for identifying patients that could potentially benefit from adjuvant chemotherapy. In addition to colon cancer, the TMRS panel was also revealed to be a reliable tool for prognostic prediction and chemotherapeutic decision-making in gastric cancer. Furthermore, we also found that the TMRS panel enabled prediction of anti-PD-L1 and anti-PD-1 immunotherapy outcomes in urothelial carcinoma patients and melanoma patients.

Implications of all the available evidence

The TMRS gene panel represents a potentially robust tool for survival prediction and treatment guidance in patients with stage I–III colon cancer and may also be applicable to other types of cancers.

Alt-text: Unlabelled Box

Introduction

Colon cancer is a major source of morbidity and mortality worldwide [1]. Currently, the AJCC staging system and histologic classification remain the most important guidelines for stratifying patients and making clinical decisions [2]. However, due to the high levels of heterogeneity found in colon cancer, prognoses may vary widely between patients with similar clinical features. Therefore, in order to stratify patients more precisely, it becomes necessary to consider other prognostic factors in addition to clinical factors.

The tumour microenvironment (TME), consisting of multiple immune cells and stromal cells, is critical for the regulation of cancer initiation and development as well as cellular response to chemotherapy [3,4]. In recent years, the rise of immunotherapy, including immune checkpoint inhibitors, has shown assessment of TME landscape heterogeneity and reshaping of immune microenvironment to be promising avenues for future cancer management [5]. In colon cancer, assessing amounts of tumour-infiltrating lymphocytes based on IHC staining of CD3 and CD8, via an “immunoscore”, is considered an important supplemental marker in the TNM staging system for relapse and mortality prediction [6,7]. Unfortunately, the accuracy of prognosis prediction using the immunoscore system was found to be limited. A multi-central study conducted by Galon et al. [6] indicated that the c-indexes of the immunoscore for relapse-free survival (RFS) and overall survival (OS) were 0·62 and 0·58, respectively. Therefore, a novel TME-related prognostic model may be needed to improve predictive accuracy and identify patients, for whom chemo- and immunotherapies may be more beneficial.

Recent advances in high-throughput gene testing technology have provided an opportunity to define the genetic landscape of colon cancer and led to the development of many molecular signatures for prognosis prediction and personalisation of treatment paradigms [[8], [9], [10], [11]]. However, some of these signatures were frequently ill-defined, having been generated from unspecified genetic backgrounds. To the best of our knowledge, prognostic signatures based on TME-relevant genes in colon cancer have not yet been proposed. The objective of the current study was to develop a robust prognostic gene panel utilising the machine learning method, designated as the “tumour microenvironment risk score (TMRS)”, to improve the risk stratification of patients with stage I–III colon cancer. Furthermore, we assessed the ability of this panel to predict patient response to chemotherapy and immune-checkpoint inhibitors.

Materials and methods

Transcriptome data acquisition and pre-processing

The Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) was searched for eligible colon cancer datasets that fulfilled the following criteria: samples were hybridised to the Affymetrix HG-U133 Plus 2·0 (GEO accession number GPL570) platforms; >50 stage I–III colon cancer patients were included in each dataset; and information on the TNM stage was available. In order to explore the role of our model in gastric cancer, we downloaded the “GSE62254” dataset, simultaneously containing RFS and OS information generated via the GPL570 platform. Finally, for the immunotherapeutic efficiency analysis, two transcriptomic datasets from patients with metastatic urothelial cancer (mUC) treated with anti-PD-L1 agents (atezolizumab, IMvigor dataset, retrieved via R software using “IMvigor” package) [12], and patients with metastatic melanoma treated with anti-PD-1 agents (pembrolizumab or nivolumab, GSE78220, downloaded from GEO website) were downloaded. Expression profiles of these two cohorts were generated via high throughput sequencing. Raw “CEL” files of microarray data were downloaded and normalised using a robust multiarray averaging method with “affy” and “simpleaffy” packages [13]. RNA sequencing data was transformed using the “voom” algorithm in order to convert count data to values similar to those resulting from microarrays [14]. The “ComBat” algorithm was applied to reduce the likelihood of batch effects from non-biological technical biases. A summary of the information of all datasets used in this study is provided (Supplement Table S1).

Study population and clinicopathological variables

Samples were randomly separated into training and validation (7:3) sets for prognostic analyses based on cohorts, in order to identify and evaluate the models as we described before using “createDataPartition” function of the “caret” package [15]. Clinical information was retrieved using the “GEOquery” package for GEO datasets and the “IMvigor” package for the IMvigor cohort. The endpoints analysed in this study were RFS, defined as the interval between the date of diagnosis and date of tumour relapse, and OS, defined as the interval between the date of diagnosis and death.

Robust tumour microenvironment prognostic gene identification

TME-relevant genes were obtained from 12 published studies [[16], [17], [18], [19], [20], [21], [22], [23], [24], [25], [26], [27]], which provided transcriptomic signatures for multiple immune and stromal cell populations. Robust prognostic genes were identified using two steps: first, we assessed the correlation between the relative expression value (z-transformed) of each gene and RFS via Cox univariate regression analysis in the entire cohort, where genes with P < 0·05 were selected for further analysis; next, we used bootstrapping to test the genes which passed initial filtering for robustness as follows: 70% patients randomly extracted from the training cohort were assessed for survival impact of their genes. This procedure was repeated 1000 times and the genes that were incorporated in 70% of resample runs (achieved P < 0·05 in robustness testing) were considered as robust prognostic genes and selected for further analysis.

TMRS gene panel generation using LASSO Cox regression

The Cox regression model, with least absolute shrinkage and selection operator (LASSO) penalty, was implemented to reduce dimensionality and select the most useful prognostic markers among the robust prognostic genes identified before [28,29]. Notably, all gene expression values were dichotomised before entering the LASSO Cox regression, and the “surv_cutpoint” function of the “survminer” R package was used to determine the optimal cut-off point of each gene based on the maximally selected log-rank statistics. Moreover, we set the “minprop” parameter of the “surv_cutpoint” function (referring to the minimal proportion of observations per group) to 30% to avoid the occurrence of too few patients in a certain group. Genes represented by optimal values of the penalty parameter λ, which were determined by ten-fold cross-validations, constituted the TMRS panel in this study. Risk scores based on the TMRS panel were also constructed using the dichotomised expression value of selected genes (a value of one or two was given to represent an expression value higher or lower than the cut-off value) via a Cox regression analysis in the training cohort. Among these, the risk score was named TMRS-RFS when RFS was used as the endpoint variable to generate the cut-off values and analysed by Cox regression, and correspondingly, the risk score was named TMRS-OS when OS was used as the endpoint.

Estimation of immune infiltration

An immune infiltration estimation was conducted using the “Microenvironment Cell Populations-counter (MCP-counter)” method, which allows robust quantification of the absolute abundance of ten immune and stromal cell populations in heterogeneous tissues from transcriptomic data [18]. The R package “MCPcounter” was applied to achieve the transformation of mRNA data to the level of non-tumour cell infiltration in the tumour microenvironment, and gene expression profiles were prepared using standard annotation files prior to MCP-counter analysis. The stromal score, immune score, and tumour purity were estimated by applying the “ESTIMATE” R package developed by Yoshihara et al. [30](https://sourceforge.net/projects/estimateproject/).

Gene set variation analysis (GSVA)

GSVA is a gene set enrichment method that estimates the variation of pathway and biological process activity over a sample population in an unsupervised manner [31]. The gene set files of “c2.cp.kegg.v6.2.symbols” and “h.all.v6.2.symbols,” downloaded from the “Molecular Signatures Database,” were employed for GSVA using “GSVA” packages for R. The significance threshold was set at an adjusted P < 0·05.

Statistical analysis

The normality of the variables was tested via the Shapiro-Wilk normality test for comparisons of two groups. The statistical significance of differences between normally distributed variables was estimated using the unpaired Student's t-test, and non-normally distributed variables were analysed via the Mann-Whitney U test. For comparisons of more than two groups, Kruskal-Wallis and one-way ANOVA tests were used as non-parametric and parametric methods, respectively. Correlation was computed using Spearman's and distance correlation analyses. Survival rates were calculated using the Kaplan–Meier method, and the significance of differences between survival curves was determined using the log-rank test. In regard to the heterogeneity between different types of cancers, the best cut-off values for each continuous prognostic marker were recalculated using the “survminer” R package separately for different tumour types. Uni- and multivariate analyses were performed using Cox proportional hazard models with the stepwise method “LR forward”. Nomogram construction and validation were performed using Iasonos' guide [32]. Survival predictive accuracy of prognostic models was assessed based on a time-dependent receiver operating characteristic curve (ROC) analysis and Harrell's concordance index (c-index) analysis. All statistical analyses were conducted using the R software (version 3.5.0) and SPSS software (version 25.0) and P values were two-tailed. Statistical significance was set at P < 0·05.

Results

Colon cancer patient characteristics and robust prognostic gene identification

A summary of the information of all datasets used in this study is provided (Supplemental Table S1). Detailed patient characteristics are listed (Supplemental Table S2). A total of 990 patients diagnosed with stage I–III colon cancer from five GEO datasets (GSE17538, GSE33113, GSE37892, GSE38832, and GSE39582) were retrospectively analysed in this study. The median age at diagnosis was 69·0 years (range, 22·0–97·0 years) and 481 (48·6%) of the patients were male. Among them, RFS information was available for 990 patients, wherein mean survival was 146 months. OS information was documented in 678 patients and the mean survival was 130 months. Through the 2-step analysis described in the “Materials and methods”, 1746 of 5952 genes passed the first filter and 797 genes, the expression levels of which were stably and significantly correlated with prognosis, were eventually identified and defined as robust prognostic genes (Supplemental Table S3).

Construction of molecular subgroups using TME-relevantrobust prognostic genes

First, we used unsupervised clustering methods in order to classify 990 tumour samples into different molecular subgroups based on 797 robust prognostic genes. The “ConsensusClusterPlus” package was used to evaluate clustering stability and select the optimal cluster number. Two distant patient clusters, termed as the tumour microenvironment cluster 1 (TMEC1) and TMEC2, were finally identified (Fig. 1a–b), and comparison of the proportion of patients from different GEO series between two TMEC clusters showed no significant differences (Supplemental Table S4). Via the log-rank test, the Kaplan–Meier curve indicated significant survival differences between the two clusters for both OS (Fig. 1c) and RFS (Fig. 1d). We further focused on GSE39582 datasets, which provided the most comprehensive patient information, characterising biological and clinical differences among these clusters. Samples in TMEC2 exhibited a more advanced tumour stage and proficient mismatch repair (pMMR) status (Fig. 1e and Supplemental Table S5). No significant distribution difference was found in terms of the CpG island methylator phenotype, chromosome instability status, and genetic mutations (KRAS, BRAF, and P53). With regards to biological behaviour (Fig.1e), pathways involved in DNA replication and repair, oxidative phosphorylation (OXPHOS), cell cycle, and one carbon metabolism were activated in TMEC1 (favourable survival), whereas TMEC2 (worse survival), as expected, showed enrichment of pathways related to stromal activation and cancer progression, such as angiogenesis and epithelial-mesenchymal transition (EMT). As confirmed by subsequent cell infiltration analysis, TMEC2 patients showed an obvious increase in infiltration by stromal cells, including neutrophils, endothelial cells, and fibroblasts (Fig. 1f). Interestingly, it was also observed that in TMEC2 patients, the enrichment of immune relevant pathways was accompanied by an increase in the expression of immune checkpoint genes. Finally, the distribution of TMEC relative to that of other established colon cancer molecular subtypes was compared. The result demonstrated that TMEC2 patients were mainly concentrated in the C4, C6, and CMS4 subtypes, mostly representing mesenchymal phenotypes, while TMEC1 patients were mainly concentrated in the C1, C3, and CMS2 subtypes, mostly displaying epithelial phenotype characteristics (Fig. 1g, Supplemental Table S5) [33,34].

Fig. 1

Consensus clustering of tumour microenvironment (TME) genes in colon cancer. (a) Consensus matrices of colon cancer patients for k = 2; (b) Colon cancer cases are divided into two subtypes based on unsupervised analysis and hierarchical clustering of 797 robust prognostic genes. Months of relapse-free survival and relapse status (relapse, red; censor, light yellow) are indicated above the heatmap; (c–d) Differences in patient overall survival (c)and relapse-free survival (d) with two clusters; (e) Heatmap showing the activation status of the biological processes in different TME-relevant clusters; (f) Violin plot of the comparison of immune and stromal cell infiltration between the different TME-relevant clusters; (g) Sankey chart displaying the distribution of the TME-relevant clusters in C1–C6 subtypes and CMS subtypes. TMEC, tumour microenvironment cluster;MMR, mismatch repair; dMMR, deficient mismatch repair; pMMR, proficient mismatch repair; CIMP, CpG island methylator phenotype; CIN, chromosome instability; MT, mutant type; WT, wild type; IC, immune check point; EMT, epithelial-mesenchymal transition; CTL, cytotoxic lymphocyte; NK, natural killer; DC, dendritic cell; EC, endothelial cells;CMS, consensus molecular subtypes.

Construction of prognostically relevant TMRS gene panel

In order to develop a gene panel practical for clinical use, we applied the LASSO Cox regression model to the 797 robust prognostic genes for dimension reduction. Patients were randomly regrouped into training and validation cohorts for prognostic analyses as described in the “Materials and methods” section. Comparison of patient characteristics between the two groups showed no significant differences (Supplemental Table S2). Through the LASSO model (Supplemental Fig. S2), we generated a TMRS gene panel consisting of 100 genes (Supplemental Table S6) and built two prognostic models using Cox analysis based on RFS (TMRS-RFS) and OS (TMRS-OS) information separately. Patients were stratified into two groups based on TMRS-RFS and TMRS-OS values, respectively, using a cut-off value calculated in the entire cohort (2·26 for TMRS-RFS and 3·02 for TMRS-OS). In both training and validation sets, Kaplan–Meier curves indicated that patients in the high-TMRS-RFS group had a significantly higher risk of relapse (Fig. 2a–b). In ROC (Fig. 2a–b) and c-index analyses (Table 1), the TMRS-RFS model showed a much higher predictive ability than that of the TNM stage, considered as a continuous variable, in the training and validation cohorts. Since stage is a categorical variable, we converted TMRS-RFS from continuous to three-classified variables to enhance comparability, and the superior predictive accuracy of TMRS-RFS was sustained even as a categorical variable (Supplemental Fig. S3a–c and Table 1). Using multivariate analysis, the TMRS-RFS model was also found to be a strong independent risk factor when treated as a continuous variable in all patient cohorts (Table 2). Similar results were also found for the TMRS-OS model in 678 patients with documented OS information (Fig. 2c, Supplemental Fig. S3d–e, Table 1, Table 2).

Supplementary Fig. S2

Feature selection for building TMRS. Ten-time cross-validation for tuning parameter selection in the LASSO model. Solid vertical lines represent partial likelihood deviance ± standard error (SE). The dotted vertical lines are drawn at the optimal values by minimum criteria and 1-SE criteria. The partial likelihood deviance versus log (λ) was plotted, where λ is the tuning parameter.

Fig. 2

TMRS panel is a prognostic marker. (a–b) Kaplan–Meier curves (left) and ROC curves (right) of relapse-free survival according to TMRS-RFS groups in the training cohort (a) and validation cohort (b); (c) Kaplan–Meier curves (upper) and ROC curves (down) of overall survival according to TMRS-OS groups; (d–e) Forest plots of the associations between TMRS-RFS and relapse-free survival (d) and the associations between TMRS-OS and overall survival (e) in various subgroups. Unadjusted HRs (boxes) and 95% confidence intervals (horizontal lines) are depicted. TMRS, tumour microenvironment risk score; RFS, relapse-free survival; OS, overall survival; HR, hazard ratio; CI, confidence interval; CMS, consensus molecular subtypes.

Table 1

Harrell's concordance indexes of the TMRS model and stage in different cohorts.

Cohort	TMRS-RFSa	TMRS-RFSb	TMRS-OSa	TMRS-OSb	Stage6th
RFS
Training	0.86 ± 0.02	0.82 ± 0.02	/	/	0.64 ± 0.02
Validation	0.80 ± 0.04	0.76 ± 0.03	/	/	0.61 ± 0.03
OS
Entire	/	/	0.80 ± 0.02	0.75 ± 0.02	0.56 ± 0.02

Abbreviation: TMRS, tumour microenvironment risk score; RFS, relapse-free survival; OS, overall survival.

Continuous variables.

Category variables.

Supplementary Fig. S3

Predictive accuracy of the TMRS panel as category variables. (a) TMRS-RFS was triple classified based on cut-off values calculated using X-tile; (b–c) Kaplan–Meier curves (left) and ROC curves (right) of relapse-free survival according to TMRS-RFS groups in the training cohort (b) and validation cohort (c); (d) TMRS-OS was triple classified based on cut-off values calculated using X-tile; (e) Kaplan–Meier curves (left) and ROC curves (right) of overall survival according to TMRS-OS groups. TMRS, tumour microenvironment risk score; RFS, relapse-free survival; OS, overall survival.

Table 2

Univariate and multivariate survival analyses of TMRS-RFS, TMRS-OS and clinical variables.

	UVA (RFS)		UVA (OS)		MVA (RFS)				MVA (OS)
	Entire	p-Value	Entire	p-Value	Training	p-Value	Validation	p-Value	Entire	p-Value
Agea	1.00 (0.99–1.01)	0.996	1.04 (1.02–1.05)	< 0.001	NE		NE		1.03 (1.01–1.04)	< 0.001
Gender (vs. Male)	0.73 (0.59–0.90)	0.003	0.82 (0.62–1.09)	0.176	NE		NE		0.66 (0.48–0.90)	0.009
TMRS-RFSa	2.32 (2.13–2.53)	< 0.001	/		2.70 (2.33–3.13)	< 0.001	1.71 (1.43–2.04)	< 0.001	/
TMRS-OSa	/		2.73 (2.38–3.12)	< 0.001	/		/		2.67 (2.30–3.10)	< 0.001
Stage(vs. stage I)						0.015	NE		NE
Stage II	8.24 (2.03–33.39)	0.003	1.74 (0.91–3.34)	0.097	5.27 (0.72–38.64)	0.101
Stage III	17.97 (4.45–72.55)	< 0.001	2.30 (1.20–4.41)	0.012	7.00 (0.96–50.97)	0.055
CMS (vs. CMS4)					NE		NE		NE
CMS1	0.55 (0.35–0.86)	0.009	0.96 (0.63–1.45)	0.846
CMS2	0.59 (0.41–0.84)	0.003	0.63 (0.44–0.92)	0.015
CMS3	0.422 (0.25–0.73)	0.002	0.38 (0.20–0.72)	0.003

Abbreviation: TMRS, tumour microenvironment risk score, UVA, univariate analysis, MVA, multivariate analysis; RFS, relapse-free survival; OS, overall survival; CMS, consensus molecular subtypes; NE, not enter.

Continuous variable.

The predictive power of the TMRS-RFS and TMRS-OS models was next tested in various subgroups stratified by patient dataset, age, gender, stage, tumour site, and CMS subtype in the entire cohort, respectively, where TMRS-RFS and TMRS-OS were both analysed as continuous variables. Forest plots indicated that, for both models, a higher value could significantly identify the patients with worse prognoses in all subgroups (Fig. 2d–e).

TMRS panel predicts therapeutic benefit of chemotherapy in colon cancer

Adjuvant chemotherapy (ADJC) is the main treatment strategy for non-metastatic colon cancer patients [35]. Since only the GSE39582 dataset recorded chemotherapy information of patients, we analysed the relationship between the TMRS panel and ADJC benefits in this dataset, where OS was used as the treatment outcome. Survival benefits of low-TMRS-RFS and low-TMRS-OS were maintained regardless of ADJC conduction (Fig. 3a–b). Interestingly, it was observed that ADJC significantly reduced the mortality risk of patients only in low-TMRS-RFS and low-TMRS-OS groups but did not confer survival benefits to patients in high-TMRS-RFS or high-TMRS-OS groups (Fig. 3c). Furthermore, the results of stratified analysis of each stage (Supplemental Fig. S4), showed that treatment benefits of ADJC were higher for patients in groups with low scores, in either Stage II or III. To develop a clinically relevant quantitative method for predicting the probability of patient mortality, we constructed two nomograms (Fig. 3d–e) integrating both TMRS panel derived scores and independent clinical prognostic factors in the GSE39582 dataset (Supplemental Table S7). The calibration plots showed that the derived nomograms performed well compared to the performance of an ideal model (Fig. 3f–g). Decision curve analysis revealed that clinical usefulness of the nomograms significantly overwhelmed the TNM stage (Fig. 3h–i).

Fig. 3

(a–b) Kaplan–Meier curves of overall survival for patients in subgroups stratified by both TMRS-RFS (a)/TMRS-OS (b), and receipt of adjuvant chemotherapy; (c) Forest plot showing the benefit of adjuvant chemotherapy in different TMRS groups; (d–e) Nomograms for predicting the probability of patient mortality based on TMRS-RFS (d), TMRS-OS (e) and clinical variables; (f–g) Plots depict the calibration of nomograms based on TMRS-RFS (f) and TMRS-OS (g) in terms of agreement between predicted and observed 2-year, 3-year, and 5-year outcomes. Nomogram performance is shown by the plot, relative to the 45-degree line, which represents the ideal prediction; (h–i) Decision curve analyses of the nomograms based on TMRS-RFS (h) and TMRS-OS (i) for 2-year, 3-year, and 5-year risk. ADJC, adjuvant chemotherapy; TMRS, tumour microenvironment risk score; RFS, relapse-free survival; OS, overall survival; Pr, probability; Nomo, nomogram.

Supplementary Fig. S4

Forest plot showing the benefit of adjuvant chemotherapy in different TMRS groups stratified by each stage. TMRS, tumour microenvironment risk score; RFS, relapse-free survival; OS, overall survival; CI, confidence interval.

Identification of TMRS-RFS and TMRS-OS related biological pathways and processes

The correlations between TMRS panel derived scores with clinical characteristics and molecular subtypes were further investigated in the GSE39582 series (Fig. 4a–b). In terms of clinical characteristics, both TMRS-RFS and TMRS-OS were increased in patients with more advanced stages and patients who had relapsed and died due to the disease. Furthermore, while gender influenced the value of TMRS-RFS, that of TMRS-OS varied between age and tumour site. In terms of molecular characteristics, we observed that KRAS mutation simultaneously up-regulated the values of TMRS-OS and TMRS-RFS, and patients in molecular subtypes C4, C6, and CMS4 exhibited significantly higher values of TMRS models than others. However, mismatch repair status was significantly correlated only with the TMRS-RFS level. Next, we used GSVA to study the association between potential biological phenotypes and the two TMRS models. The result showed that the biological behaviour of TMRS-OS and TMRS-RFS was similar to each other to some extent (Fig. 4c). The two models were both positively correlated with stromal activation relevant pathways, but negatively correlated with those of DNA replication and repair, OXPHOS, and one carbon metabolism. Cell infiltration analysis indicated that infiltration by multiple stromal cell types significantly increased with rising TMRS-OS and TMRS-RFS scores, while some immune cells, such as cytotoxic lymphocytes and NK cells, were negatively correlated with TMRS-RFS scores (Fig. 4d).

Fig. 4

Clinical significance and biological function of TMRS panel. (a–b) TMRS-RFS (a) and TMRS-OS (b) values in different clinical subgroups. Boxes represent 25–75% of values, black lines in boxes represent median values, whiskers represent 1.5 interquartile ranges, and black dots represent outliers; (c) Correlation matrix of TMRS-RFS, TMRS-OS values and the activation levels of biological process. Shading colour represents the value of corresponding correlation coefficients and nonsignificant correlations are denoted by “x”; (d) Bubble plot showing the correlations between TMRS-RFS, TMRS-OS values and cell infiltration. TMRS, tumour microenvironment risk score; RFS, relapse-free survival; OS, overall survival; RS, relapse status; LS, live status;MMR, mismatch repair; CIMP, CpG island methylator phenotype; CIN, chromosome instability; MT, mutant type; WT, wild type; CS, c subtype; CMS, consensus molecular subtypes; APP, antigen processing and presentation; HCL, hematopoietic cell lineage; TCR, T cell receptor; TCA, tricarboxylic acid; PM, pyruvate metabolism; OXPHOS, oxidative phosphorylation; PUM, purine metabolism; PYM, pyrimidine metabolism; CYS, cysteine; OC, one carbon; BER, base excision repair; NER, nucleotide excision repair; MMR, mismatch repair; HR, homologous recombination; NHEJ, non-homologous end joining; EMT, epithelial-mesenchymal transition; NK, natural killer; NS, not significant.

The TMRS panel could be used for prognostic prediction and tailoring therapies in gastric cancer

Reportedly, TME plays an important role in the initiation and progression of multiple tumour types [3,4]. Therefore, the current study investigated whether the TMRS panel can predict prognoses in patients with stage I–III gastric cancer. The GSE62254 series downloaded from GEO datasets indicated that in both OS (Fig. 5a) and RFS (Fig. 5b), TMRS panel derived scores were predictive of a poor outcome and acted as independent prognostic factors in gastric cancer patients (Supplemental Table S8). A subsequent ROC analysis further confirmed that the TMRS panel derived score models exhibited much better predictive accuracy in gastric cancer compared to those of the TNM stage (Fig. 5a–b). Additionally, we found that the TMRS-RFS model, but not the TMRS-OS model, could be used to screen gastric cancer patients who could benefit from ADJC (Fig. 5c–d). Subsequent cell infiltration (Fig. 5e) and GSVA analyses (Fig. 5f) showed that, compared with TMRS-OS model, the TMRS-RFS model in gastric cancer is markedly more related to stromal cell infiltration and the activation level of stromal-relevant biological processes, such as EMT and angiogenesis. Finally, the distribution of TMRS scores in patients of different ACRG subtypes was compared [36]. The result showed that both TMRS-OS and TMRS-RFS models have the highest scores in the “EMT” subtype, while the lowest scored were observed in the “MSI” subtype (Fig. 5g).

Fig. 5

Exploration of the role of the TMRS panel in gastric cancer cohort. (a–b) Kaplan–Meier curves (left) and ROC curves (right) of overall survival according to TMRS-OS groups (a) and relapse-free survival according to TMRS-RFS groups (b); (c–d) Kaplan–Meier curves of overall survival for gastric cancer patients in subgroups stratified by both TMRS-OS (c), TMRS-RFS (d), and receipt of adjuvant chemotherapy; (e) Bubble plot showing the correlations between TMRS-RFS, TMRS-OS values and cell infiltration in GSE62254 series; (f) Heatmap showing the correlation between TMRS-RFS, TMRS-OS and biological process in GSE62254 series; (g) Violin plot displayed the distribution of the TMRS-OS value (up) and TMRS-RFS value (down) in different ACRG molecular subtypes of gastric cancer. P value for TMRS-OS comparison marked in red, and TMRS-RFS comparison marked in yellow. Boxes represent 25–75% of values, black lines in boxes represent median values, white dots represent mean values, and whiskers represent 1.5 interquartile ranges. TMRS, tumour microenvironment risk score; RFS, relapse-free survival; OS, overall survival; ADJC, adjuvant chemotherapy; NK, natural killer; NS, not significant; OXPHOS, oxidative phosphorylation; TCA, tricarboxylic acid; APP, antigen processing and presentation; TCR, T cell receptor; BER, base excision repair; NER, nucleotide excision repair; MMR, mismatch repair; HR, homologous recombination; NHEJ, non-homologous end joining; EMT, epithelial-mesenchymal transition; MSI, microsatellite instability; MSS, microsatellite stability. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

TMRS panel predicts immunotherapeutic benefit

Currently, effective predictive markers for immune therapy are limited. The identification of novel predictive markers is crucial for further advancing precision immunotherapy. We retrieved two immunotherapy datasets (Imvigor210, GSE78220) with the transcriptome data released to explore whether TMRS panel could predict immunotherapeutic benefit. Imvigor210 documented expression data in human mUC samples from patients who did or did not respond to anti-PD-L1 immunotherapy, and GSE78220 is a malignant melanoma dataset of patients who received anti-PD-1 therapy. The Kaplan–Meier curve (Fig. 6a) revealed that a higher TMRS-OS value was associated with much poorer survival in mUC patients who received immunotherapy. A violin plot further demonstrated that TMRS-OS values successively and significantly increased in patient groups with complete or partial response to stable disease and progressive disease(Fig. 6b). In total, 57 (46·3%) responders of immunotherapy in the low-risk group and 11 (6·3%) responders in the high-risk group (P = 3·0 × e−16, Chi square test, Fig. 6c) were identified. Notably, the TMRS-OS model enabled the classification of responders and non-responders with a considerably higher predictive power (AUC = 0·84, Fig. 6d) than both neoantigen (AUC = 0·77) and tumour mutation burden (TMB, AUC = 0·73). A complex model, combining the TMRS-OS, neoantigen, and TMB models using logistic regression, identified treatment response with a level of accuracy as high as 88% (Fig. 6d). PD-L1 expression in immune cells, but not tumour cells, was associated with the response of mUC patients towards anti-PD-L1 therapy [12]. Coincidently, a significant negative correlation between the TMRS-OS value and PD-L1 expression in immune cells was also observed (Fig. 6e). In addition, immune phenotype analysis demonstrated that the TMRS-OS value decreased with the activation of inflammation, although the P value was not significant (P = 0·280, one-way ANOVA tests). The correlation test showed a negative correlation between TMRS-OS and immune cell infiltration, TMB, and neoantigen, respectively. Unfortunately, infiltration of stromal cells was not significantly associated with TMRS-OS value (data not shown). Finally, the GSE78220 dataset was used to further verify predictive effect of TMRS panel on treatment outcome of melanoma immunotherapy. The result also showed that the TMRS panel could discriminate the prognoses of different patients (Supplemental Fig. S5a), and that patient scores were well connected to the treatment effect (Supplemental Fig. S5b–c), with the AUC reaching a value of 0·86 (Supplemental Fig. S5d). Of note, since information of only 27 patients was included in the GSE78220 series, the use of the TMRS panel containing 100 genes for the scoring might result in an overfitting problem. Thus, caution must be exercised when interpreting the relevant results.

Fig. 6

TMRS panel could be used for predicting immunotherapeutic benefit. (a) Kaplan–Meier curves of overall survival according to TMRS-OS groups in IMvigor210 cohort; (b) Violin plot showing the correlation of TMRS-OS value and immunotherapeutic effectiveness in the Imvigor cohort. Boxes represent 25–75% of values, black lines in boxes represent median values, white dots represent mean values, and whiskers represent 1.5 interquartile ranges; (c) Rate of clinical response to anti-PD-L1 immunotherapy in high or low TMRS-OS groups. (d) ROC curves of immunotherapy response predictions according to different biomarkers; (e) Box plot displaying the distribution of TMRS-OS values in patients with different immune cell and tumour cell PD-L1 status, as well as immune phenotypes. Boxes represent 25–75% of values, black lines in boxes represent median values, whiskers represent 1.5 interquartile ranges, and black dots represent outliers; (f)The correlation chord chart showing the mutual correlation between TMRS-OS value, TMB, NEO, and immune cell infiltration. TMRS, tumour microenvironment risk score; OS, overall survival; CR, complete response; PR, partial response; SD, stable disease; PD, progressive disease; NEO, neoantigen, TMB, tumour mutation burden; IC, immune cell; TC, tumour cell; CTL, cytotoxic lymphocyte; NK, natural killer.

Supplementary Fig. S5

The role of TMRS panel in predicting immunotherapeutic benefit in GSE78220 series. (a) Kaplan–Meier curves of overall survival according to TMRS-OS groups; (b)Violin plot showing the correlation of TMRS-OS value and immunotherapeutic effectiveness. Boxes represent 25–75% of values, black lines in boxes represent median values, white dots represent mean values, and whiskers represent 1.5 interquartile ranges; (c) Rate of clinical response to anti-PD-1 immunotherapy in high or low TMRS-OS groups; (d) ROC curves of immunotherapy response predictions according to TMRS-OS. TMRS, tumour microenvironment risk score; OS, overall survival; CR, complete response; PR, partial response; PD, progressive disease.

Discussion

The heterogeneity of TME contains multiple dimensions of information on patient prognosis and treatment response. Currently, several signatures based on TME have been reported [20,21,37]. However, as the genes included in these signatures were not all correlated with prognosis, the prognostic prediction ability of these signatures were not satisfactory. In colon cancer, Galon et al [38]. developed an immunoscore system based on the density of CD3+, CD8+, or CD45RO+ lymphocytes in the central- and peri- tumoral areas, as represented by the intensity of IHC staining, and found that the prognostic ability of the immunoscore was stronger than, and could potentially replace, the TNM stage. However, the latest multi-central clinical research reported that the accuracy of the immunoscore was only approximately 60% [6]. As such, although the prognostic ability of the immunoscore is higher than that of the TNM stage, there is still room for improvement in prediction accuracy. Neglecting the assessment of stromal cell infiltration may limit immunoscore prediction accuracy. When comparing the cell infiltration status of TMEC1 and TMEC2, we found that both stromal cell and immune cell infiltration significantly increased in the TMEC2 type, and the level of increase in the former was far higher than the latter. Further, the activation level of some stromal pathways and expression of immune checkpoint molecules were also significantly higher in the TMEC2 type. This indicated that large amounts of stromal cell infiltration not only reversed the survival benefit of immune infiltration, but also kept the infiltrated immune cells in a state of exhaustion, with inhibited functions, through up-regulation of immune checkpoint molecules [39]. Under such circumstances, the prediction accuracy of immunoscore may greatly decrease. In contrast to the immunoscore, the TMRS panel we developed was based on genes from 12 previously published studies, which provided transcriptomic signatures of multiple non-cancerous cell types located in TME. These cells included not only immune cells of various functional subtypes, but also multiple stromal cells, such as endothelial cells, adipose cells, and fibroblasts. Meanwhile, the successive application of the bootstrap method to identify robust prognostic genes, and the machine-learning method, LASSO regression, to screen the optimal combination of genes, markedly raised the accuracy of the TMRS panel in predicting relapse and mortality risks in colon cancer patients. Thus, this model may have a strong clinical transformation value. However, microarrays based on whole transcriptomes are not practical clinically. According to the guidelines established by Altman et al. [40], before the gene signature is applied as a clinical grade assay, identification of an appropriate approach to quantify the expression (microarray) and the design of specific probes based on the sequences tested in the microarray chips are two necessary steps. Therefore, to make the TMRS signature more clinically applicable, we plan to use a nanostring technique to develop a gene detection kit based on these 100 genes and re-standardise the gene expression value and the cut-off value. Finally, only signatures validated in independent cohorts of patients with full clinical annotation could be applied clinically, which will be reported in a future paper. We will first validate the prognostic value of this model at our centre.

Chemotherapy is currently the main treatment strategy for cancer. Detecting patients that could potentially benefit from chemotherapy is an important step in precision treatment. Through analysis of the GSE39582 dataset, we found that the TMRS panel could effectively screen patients who were responsive to ADJC. Patients in the lower TMRS score group may benefit significantly more than patients in the higher TMRS score group. Several studies have reported that the EMT and angiogenesis are the main factors influencing the efficacy of chemotherapy [[41], [42], [43]]. Correspondingly, we found that the above-mentioned pathways were significantly more activated in colon cancer patients with higher TMRS scores than in patients with lower scores. This finding has also been confirmed in gastric cancer datasets where TMRS-RFS, with significant relevance to stromal pathway activation level, could be used to discriminate patients who could benefit from ADJC, while TMRS-OS, with no such significant relevance, did not display this function. Of note, due to limitation of the sample size in the datasets, additional prospective studies are required to further verify these findings.

Recently, clinical results from trials investigating checkpoint blockade inhibitors have attracted high interests in this therapeutic modality [5]. However, the selection criteria for candidates who are likely to benefit from such regimens requires further investigation. Encouragingly, the risk score generated by TMRS panel was also found to be predictive of outcome after anti-PDL1 in mUC and anti-PD1 in melanoma. The prediction accuracy of both patient cohorts reach approximately 85%, which is significantly higher than the tumour mutation burden and abundance of neoantigens models, which are biomarkers that predict patient's response to immunotherapy [44,45]. Currently, in colon cancer patients, the mismatch repair status is the only marker which predicts whether the patient should receive immunotherapy, and only patients with dMMR could possibly benefit [46]. As our study also revealed significant variation of the TMRS-RFS value between patients with different mismatch repair status, it may be speculated that immunotherapy might also be a preferable choice for patients in the low-TMRS-RFS group. Further investigation on the association between TMRS panel and immunotherapy efficiency in colon cancer patients are warranted.

There were some limitations to the present study. Firstly, the patient population was heterogeneous because of the retrospective nature of this study. Secondly, all colon cancer transcriptome profiling used for panel construction was produced by the GPL570 platform. Therefore, caution should be exerted when applying the panel to samples tested using platforms other than GPL570. Thirdly, since gene expression data was entered into Cox regression as categorical variables, the optimal cut-off value needs to be further verified in future studies.

In conclusion, the TMRS gene panel is a robust tool for survival prediction and treatment guidance in patients with stage I–III colon cancer. Further, prospective clinical trials are warranted to validate our findings.

The following are the supplementary data related to this article.

Supplemental Table S1

Basic information of gene expression profiling series used in this study.

Supplemental Table S2

Patients' basic characteristics.

Supplemental Table S3

Robust prognostic genes.

Supplemental Table S4

Comparison of the proportion of patients from different GEO series between two clusters.

Supplemental Table S5

Clinical characteristics of TMEC in GSE39582 series.

Supplemental Table S6

Genes of TMRS panel.

Supplemental Table S7

Multivariate survival analyses of TMRS-RFS, TMRS-OS and clinical variables in the GSE39582 series.

Supplemental Table S8

Multivariate survival analyses of TMRS-RFS, TMRS-OS and clinical variables in the GSE62254 series.

Workflow chart for data collection and analysis.

Funding sources

This work was supported by the National Natural Science Foundation of China (No. 81772580 to Wangjun Liao) and Guangzhou Planed Project of Science and Technology (No. 201803010070 to Wangjun Liao).

Declaration of interests

The authors declare that they have no competing interests.

Author contributions

RZ, DZ and WL contributed to the planning of the study. RZ, DZ, NH and WL drafted the manuscript and revised the manuscript. NH verified the numerical results by an independent implementation. RZ, DZ, JZ, HS, and JW prepared all the figures and tables. LL, NL, MS, JB, and YL contributed to interpretation of data and review of the manuscript. All the authors reviewed and approved the final manuscript.

Availability of data and materials

GSE17538 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE17538); GSE33113 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE33113); GSE37892 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE37892); GSE38832 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE38832); GSE39582 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE39582); GSE62254 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE62254); GSE78220 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE78220); IMvigor210 R package “IMvigor210”.