A robust classifier for the intrinsic consensus molecular subtypes in colorectal cancer.

Background
Colorectal cancer presents a complex, heterogeneous bulk gene expression profile that converges on four major molecular subtypes, previously described as CMS1 to CMS4 [1]. The stroma-rich CMS4 subtype is typically associated with worse relapse-free survival, while the CMS1 subtype is enriched in tumors with microsatellite instability [1]. These subtypes arise from a mixture of multiple cell types with very different gene expression patterns.
The advent of single-cell RNA sequencing enabled the analysis of individual cell types. An analysis focused on epithelial cells from single-cell RNA sequencing of five colorectal cancer cohorts with 373’058 cells showed two main subtypes which are strongly correlated with bulk CMS2/3 and were accordingly labelled intrinsic CMS2/3 (iCMS) [2]. The bulk CMS1 subtype is almost exclusively associated with epithelial iCMS3, while the stroma-rich bulk CMS4 subtype can be associated with either epithelial iCMS2 or iCMS3. In addition, MSI tumors are also almost always iCMS3, with obvious treatment consequences. Overall, the iCMS3 seems to follow the serrated pathway biology (gastric metaplasia) with common BRAF mutations, while the iCMS2 is more closely related with the adenoma pathway biology.
Interestingly, the epithelial component can have a profound effect on prognosis. Bulk CMS4 with iCMS3 epithelial cells has worse survival than CMS4 with iCMS2 epithelial cells, showing how the iCMS is orthogonal to stromal components. The iCMS classification of cell lines also showed differential sensitivity to numerous drugs, notably iCMS2 cell line seemed more sensitive FOLFOX, while iCMS3 cell lines more sensitive to FOLFIRI, for example [2]. In the future, the iCMS could hopefully reveal clinically meaningful differences in treatment response, as already demonstrated for CMS. [3] Finally, the iCMS are generally congruent in multiregional sequencing [4], an important feature for clinical applications where tumor heterogeneity may confound tumor classification.
Wide adoption of the iCMS classification and further research on its clinical utility would be greatly facilitated by a standardized algorithm implemented in an easily available software package. Clinical application of gene expression classifiers must be robust to technical variation and batch effect, as patients are classified individually, for example during clinical trial screening, and are not always part of a batch cohort. The previously published method for iCMS classification was not specifically optimized for application to individual samples, when batch correction and normalization is not feasible [2]. Here we provide the implementation of a single-sample classifier that replicates the results of the reference algorithm with improved performance in a new, unseen, non-batch-corrected validation dataset.

Methods

Population
The study population has been previously described [5]. Briefly, the Geneva Cancer Registry database collected data and tissue for all colorectal cancer patients diagnosed between 1985 and 2013. Samples with available fresh frozen material were selected for gene expression profiling with RNAseq (N = 384), these ranged from 2002 until 2011.

RNA sequencing
Cryosections of the fresh frozen blocks were cut and the first and last sections stained using standard techniques and the tumor content verified by a trained pathologist. Only samples with a tumor cell content of more than 30% were used for further analysis. RNA was extracted from the intervening sections using the RNAeasy kit (Qiagen, Netherlands) and its’ quality and quantity controlled using a tape station (Agilent, Santa Clara, CA, USA).
RNA-seq libraries were prepared using the TruSeq Stranded Total RNA Library Prep Human/Mouse/Rat reagents (Illumina, San Diego, USA). Following quality assessment on the Fragment Analyzer (Agilent, Santa Clara, USA), libraries were 100nt paired-end sequenced to a target of 50 million clusters per library on the Illumina NovaSeq 6000 (Illumina, San Diego, USA).
The GTR cohort has a few outlier samples, based on both RNAseq quality metrics and PCA (Figure S1A). Use of the RUV algorithm for batch effect correction seems to attenuate the difference of these samples from the main population and also slightly reduces total variance (Figure S1B) [6]. We did not remove samples that were of potentially lower quality, given our goal to develop a classifier that will be useful in a real-life setting, where sample quality can vary, and the acquisition of a new sample can incur delays and risk for patients.

Single-sample classifier
For the single-sample classifier (SSC) we focused on a subset of the previously published iCMS genes [2] whose expression in epithelial cells from public single-cell datasets [7] was statistically significantly higher than in all other major cell types by a log-fold change of at least 1 and whose log-normalized average expression was also at least 1. These genes are conceptually less likely to be sensitive to biological noise from other cell types. The final gene subset contains 201 genes, shown in Table S1.
The SSC is based on non-parametric correlation distance (Kendall’s τ) from synthetic nearest neighbors. Kendall’s τ was chosen as default because it is non parametric and has greater statistical efficiency than Spearman’s ρ, especially for small datasets [8, 9]. Alternatively, the SSC also allows the use of Pearson’s, Spearman’s and cosine distance in the function parameters.
We used three non-batch corrected datasets from TCGA [10], PETACC-3 [11] (E-MTAB-990) and GSE39582 [12] with a total of 1779 samples to generate synthetic centroids for each class (Fig. 1A). Importantly, these datasets were included during the initial iCMS derivation, and the reference iCMS class of each sample is known. The expression of the 201 genes was corrected by mean-centering for the purposes of visualization, since the three datasets are based on different assays (array and RNAseq).
Each training dataset was used to extract 32 exemplar synthetic centroids for each class, for 192 centroids in total, by calculating the mean expression of each gene in 8 randomly selected samples of the same class (iCMS2 or iCMS3). The most likely class for each input sample is estimated by comparing the 90th percentiles of the distance distributions for each class (distribution quantile [DQ] criterion) for each class with an absolute cutoff that can be set to zero if the nearest, most likely class is desired or higher if a confident call is necessary. During training, a minimum correlation distance of 0.1 or a minimum difference in correlation distance of 0.05 was used as a criterion to determine confident calls (Fig. 1B).
Alternatively, a K-nearest neighbor voting criterion (KNN) can also be used, with a default K = 10 for similar results. In this case, a confident call is made if at least 90% of neighbors (default = 9) are of the same class. A function wrapper is provided to enable parallel computation of multiple samples.
No other normalization or batch correction was applied, and significant batch effect remained, a fact that was intended to capture real life conditions and ensure that the classifier would be robust to such variation (Fig. 1C).
The previously published bulk RNAseq classification was based on nearest template prediction and the reference iCMS genes [2, 13]. A fast parallel implementation of this algorithm was developed for the purposes of performance comparisons and is provided in the accompanying R package (bulk NTP classifier).

Generating reference calls in the validation dataset
Ground truth in the unseen data (GTR cohort) was established by integrating the GTR cohort with the full previously published bulk dataset, that was used during the derivation of the iCMS [2]. Briefly, the reference dataset was clustered using hierarchical Ward’s method, as previously. Based on the ordering in the tree, each new unseen sample was assigned an integer index indicating the position and its proximity to the nearest reference sample based on Pearson’s correlation. The indices of the respective positions of the nearest references were assigned to the test samples, and used to sort the samples on the test heatmap. The genes of the test heatmap used the same subset and ordering as in the reference heatmap. This method is time and compute intensive since it uses a high number of samples (N = 2873) and genes (N ~ 12000) to establish the best matching iCMS.

Additional validation
Additional validation was performed with a public RNAseq dataset from ArrayExpress, accession E-MTAB-12862 [14]. We downloaded raw count data from 1183 samples, of which 1062 were tumor samples. Clinical data, including overall survival, histology, stage, tumor site and tumor purity were also obtained from ArrayExpress. Again, raw data were minimally processed with log transformation and no attempt at batch correction or normalization.

Gene enrichment analysis
Gene set enrichment analysis was performed using the clusterProfiler R package (v4.8.3). Hallmark gene sets were obtained from MSigDB via the msigdbr package and formatted into a two-column TERM2GENE data frame (pathway name and Entrez gene ID). The previously published iCMS genes and the selected SSC genes were converted to Entrez IDs and input as separate clusters to compareCluster(fun=“enricher”, TERM2GENE = TERM2GENE, … ), with Benjamini–Hochberg adjustment for multiple comparisons. Results were summarized and visualized using dot plots in enrichplot.

Gene Set variation analysis (GSVA)
We used the popular GSVA method to generate estimates of gene set enrichment for the iCMS2 and iCMS3 genes that we also used for the SSC. We classified every sample as iCMS2 if the estimate for iCMS2 was greater than the respective estimate for iCMS3.

Statistical methods
The analyses were performed in the R environment for statistics. Survival analyses were performed with the survival packages using the logrank test. The alpha cutoff for statistical significance was set to 0.05.

Software availability
The method presented here is available in an R package that can be downloaded from github (https://github.com/CRCrepository/iCMS.SSC).

Results

Gene selection and training
The iCMS genes were differentially expressed between two distinct epithelial subgroups, in the analysis of single-cell RNAseq colorectal cancer data [2]. The iCMS2 genes seem to align with the tubular adenoma pathway (crypt bottom), while the iCMS3 genes are more closely related with the serrated sessile lesion pathway and gastric metaplasia (crypt top). As discussed above, the genes used for the classification were selected from the iCMS genes based on strong and specific epithelial expression. In the training set, use of the epithelial genes (N = 201) demonstrated a significant increase in the desirable distance variance explained by class (Fig. 2A), while reducing unrelated variance (Fig. 2B), compared with all differentially expressed genes between iCMS2 and iCMS3 (N = 431). In addition, gene set enrichment analysis of the C2 curated signatures, C4 cancer signatures, C8 cell type signatures from MSigDB and the Reactome signatures (Fig. 2C–F) showed that the epithelial gene subset was enriched for the same pathways as the larger iCMS gene set, with favorable gene ratios.
During training, the SSC confident calls achieved over 98% agreement with the reference calls with both the DQ and KNN criteria (Fig. 2G). A confident call could be made for the vast majority of training samples: 89.2% of samples with the DQ criteria, 87.3% of samples with the KNN criteria. Agreement was also high when using the nearest (not necessarily confident) calls, at 92.2% for DQ and 93.8% for KNN. Of note, the popular GSVA method using nearest calls reached an agreement of 88.6% in the training data (Figure S1D) with the identical iCMS2 and iCMS3 gene sets, and was not further developed. As mentioned in the methods section, different distance metrics are available in the software, with Kendall’s nonparametric distance chosen as default. All metrics performed similarly in training data, but agreement with Kendall’s was numerically slightly higher (accuracy Kendall’s 93%, Spearman’s 92.6%, Pearson’s 92.9%, cosine 93%, Figure S1E).
Using an absolute difference of less than 0.1 between the closest iCMS2 and iCMS3 neighbors as cutoff led to one fourth (25.5%) of the samples bearing similarity to both classes. Although a “mixed” iCMS is not formally defined at the biological level, these samples did seem to have slightly worse RFS than the “pure” iCMS classes (at cutoff 0.1, HR = 1.25, 95% 1.05–1.48, p = 0.011, Fig. 2H). Agreement of confident calls with reference reached 100% at that cutoff, and the proportions of iCMS2 and iCMS3 remained constant across a wide range of cutoff values (Figure S1F). Given the lack of a formal “mixed” class, the default cutoff for confident calls was specified at 0.05, to maximize the number of confident calls while maintaining a very high agreement with reference (Figure S1F).

The geneva tumor registry colorectal cancer cohort
The gene expression data from the Geneva Tumor Registry colorectal cancer cohort (GTR) presented a typical colorectal cancer profile. As shown in Fig. 3A, right-sided tumors were more often MSI and CMS1, and almost all BRAF-mutated tumors were right-sided. The CMS1 group was infiltrated with immune cells (orange genes). The CMS4 group was stroma-rich (green genes), and was split into an iCMS2 and iCMS3 subgroup, as previously demonstrated.
The GTR cohort was not used during the development of the SSC but was solely used for validation. By integrating and re-clustering the data with the previously published reference dataset (N = 2873 samples, Fig. 3Ein [2]) we were able to estimate the most likely iCMS class for each GTR sample, if it had been included in the initial iCMS derivation. This class was used as the reference iCMS. Using this time-consuming and resource intensive approach, 87.6% of samples could be confidently mapped to a most probable iCMS class.
The SSC achieved excellent agreement with reference calls with both the DQ and KNN criteria, surpassing 90% agreement for confident calls (Fig. 3B). The label swaps of iCMS3 to iCMS2 were more common than the opposite, especially with the DQ criteria, while most of the ambiguity concerned samples that could not be confidently classified by either the reference or the SSC (Fig. 3C, Figure S2A-D). The fraction of calls that were considered confident was quite high for all bulk methods, over 85% without applying the RUV algorithm for batch correction (Fig. 3D).
As discussed, the SSC can always deliver a nearest sample call, even when not confident. When the most likely, not necessarily confident, calls were considered, the SSC achieved a concordance ranging from 86 to 88% with or without batch correction (Fig. 3E), outperforming the previously published bulk NTP predictor, which only reached 75.4% without batch correction (89.8% with batch correction). Both the KNN and DQ criteria were robust to batch correction and achieved very comparable performance with, and without batch corrections (Figs. 3B, 3F). Interestingly, the proportion of confident calls with the KNN criterion somewhat increased with batch correction, but the agreement of confident calls slightly decreased.
Additional validation was performed with the public E-MTAB-12862 dataset (N = 1062 tumor samples). A heatmap of the data shows a similar gene expression structure, with recognizable immune (CD3), stroma (DCN, COL1A1), Wnt (LGR5) and proliferation (MKI67) gene expression clusters (Figure S2E). Again, the SSC achieved excellent agreement with reference calls with both the DQ and the KNN criteria (Figure S2F-G), reaching an accuracy of 90.6% and 90.1% respectively for the nearest sample calls, and outperforming the bulk NTP classifier (81.4%). The iCMS calls for this dataset can be found in Table S2.

Robustness to perturbation
We also tested the resistance of the SSC and the bulk NTP methods to the addition of noise in the gene expression data, simulating situations where the results are of poor quality. Specifically, we added random uniform noise in increasing ranges defined as multiples of the standard deviation of each gene in the validation GTR cohort. This resulted in a significant fraction of the total gene expression variance being random and therefore non-informative, from 10% to more than 80% on average, with obvious changes in the gene expression heatmap (Figure S3 shows examples for -/+1SD, 2SD and 3SD). All methods, including the previously published bulk NTP classifier, were quite robust to the addition of noise and maintained very stable performance for nearest iCMS calls, up to about 70% of noise in the total variance (Fig. 4A). Performance declined rapidly thereafter, but the SSC remained superior, with a slight advantage for the KNN criteria.
Ideally, the SSC should also be applied in situations in which the number of available genes may be reduced, for example if some genes have low coverage. We therefore tested what the performance would be if random iCMS genes, selected from the 201 classifier genes, were not available for analysis (Fig. 4B). Again, the performance was remarkably stable when most of the genes were removed and remained acceptable even when only 20–30 genes were available. The bulk NTP classifier is shown for reference, although it normally relies on the larger full set of iCMS genes. Its performance was on average relatively constant without batch correction, but presented much more variation depending on which genes were available for analysis.
Inversely, in some cases the available genes may be known or chosen beforehand, for example if a gene panel is designed from the ground up for the purpose of iCMS classification. We simulated this situation by progressively removing available genes, while keeping the most significantly differentially ones. Performance declines almost linearly for both the SSC DQ and KNN as the size of the gene set is reduced gradually to 19 genes (Figure S4). Ordered gene selection seems to have a positive effect on performance below 50 genes for the SSC DQ and KNN. Interestingly, ordered gene selection may improve the performance of the bulk NTP classifier, but it remains below SSC across a wide range of gene subsets (Figure S4).
Some tumors may have similarity with both iCMS2 and iCMS3, possibly because of mixed tumor components. It is not clear how the classification would be influenced by variable iCMS2/iCMS3 mixture. We therefore tested the behavior of different methods in synthetic (artificial) mixtures of iCMS2 and iCMS3 samples. Briefly, for each synthetic sample we picked three random samples from our validation cohort and calculated the mean gene expression. Thus, we generated four cohorts of synthetic samples with pure iCMS2 content, pure iCMS3 content and mixed content (1/3 or 2/3 iCMS2, the rest being iCMS3), as shown in Figure S5A. The profile of mixed samples was mostly aligned with the proportion of iCMS2/3 content, with samples containing a majority of iCMS2 aligning with pure iCMS2 and inversely samples with 2/3 iCMS3 aligning with pure iCMS3. The likeness of the sample to iCMS2 or iCMS3 as estimated by the DQ classifier mostly reflected the composition (Figure S5B), with a slightly higher sensitivity for iCMS2, showing good calibration. The bulk NTP algorithm was very sensitive to the presence of iCMS2 content. All samples with pure iCMS2 or 2/3 iCMS2 content were labeled as iCMS2. However, only 14% of samples with 2/3 iCMS3 and only 64% of samples with pure iCMS3 content were labeled as iCMS3 (Fig. 4C).
Finally, the proportion of cancer cells (sample purity) can be quite variable in tumor samples. Information about sample purity was available in the E-MTAB-12862 dataset, allowing us to examine the accuracy of each classifier across purity bins, from 10 to 80% (Fig. 4D). Both SSC DQ and KNN have remarkably stable performance in low purity samples, keeping accuracy close to 90% and strictly above 85% in all cases, while the NTP dropped below 80% when purity was below 30%.

Correlation with clinical endpoints
The SSC was applied to unseen data from the VELOUR trial to test for associations with relapse-free and overall survival. The bulk NTP classifier was prognostic with a hazard ratio (HR) of 0.778 (0.62–0.97, p = 0.025) for PFS and a HR of 0.653 for OS (0.52–0.82, p = 0.0002), both in favor of iCMS2 (Fig. 5A–B). Similarly, the SSC iCMS was prognostic with the DQ criteria for PFS (HR 0.752, 0.60–0.94, p = 0.0112) and OS (HR = 0.601, 0.48–0.75, p < 0.0001) as shown in Fig. 5C–D, but also with the KNN criteria (PFS HR 0.702, 0.56–0.88, p = 0.0026, OS HR = 0.571, 0.45–0.72, p < 0.0001), in Fig. 5E–F. Overall, the SSC seems to perform at least as well as the NTP classifier.
The SSC was applied to the E-MTAB-12862 dataset for additional validation (Fig. 6A–C). In metastatic patients, the bulk NTP classifier was not prognostic for overall survival (HR = 0.829, 0.54–1.27, p = 0.39), while both the DQ (HR 0.569, 0.38–0.86, p = 0.007), and KNN were prognostic (HR 0.542, 0.36–0.82, p = 0.0038). The hazard ratios that we observed in this dataset were highly comparable with those that we found in the VELOUR trial. Importantly, when the full cohort was considered (N = 1062), the DQ and KNN were prognostic in multivariable models including known predictors of prognosis (stage, age, grade) and CMS, while the bulk NTP was not prognostic in this setting (Fig. 6D–F). The same was observed in a bivariable model of MSI status and iCMS (Figure S6).

Compute time requirements
The bulk NTP algorithm is permutation-based and requires a relatively long computation time. The total CPU time required to execute the bulk NTP classifier on the GTR cohort was 124.5 seconds, against 12.5 seconds for the SSC(DQ) and 12.4 for SSC(KNN). The SSC and the bulk NTP algorithm were both implemented for parallel execution. With 8 parallel jobs the real elapsed time (“wall clock” time) was shorter, at 19.8, 2.0 and 1.7 seconds respectively in a mainstream personal computer (Figure S7).

Discussion
We present a robust iCMS classifier (SSC) that is specifically tailored for analyses in a per sample setting, where batch correction methods such as RUV are not applicable. The SSC obtains results that are highly concordant with the iCMS classification when replicated by clustering new, previously unseen validation cohorts with the reference iCMS dataset. Further, it is robust to technical and biological perturbation and remains highly prognostic in a retrospective analysis of the VELOUR clinical trial and a large external cohort.
The SSC was developed with minimal processing and quality control of the training datasets, a process that resulted in significant variation between the centroids that are used to infer similarity with iCMS2 and iCMS3. Instead of using two templates as in the NTP method, the SSC extends this to 192 centroids, which should better capture technical and biological variability. Further, the reduction of the reference genes to those that are predominantly and specifically expressed in epithelial cells seems to also increase the desirable variability between classes, without affecting the representation of relevant pathways.
As shown in the initial iCMS publication, the iCMS prognostic value derives mostly from splitting the stroma-rich CMS4 samples by epithelial iCMS class, with iCMS3/CMS4 having the worst prognosis of all molecular subtype combinations. Here we show that the prognostic value of the SSC is comparable with that of the bulk NTP classifier in clinical trial data. Based on the analysis of synthetic “mixed” samples, the NTP classifier seems particularly sensitive to the presence of an iCMS2 component. It is therefore possible that for samples that are biologically heterogeneous the NTP will tend to classify them as iCMS2, a hypothesis that is compatible with the greater proportion of iCMS2 in the VELOUR trial. Further analyses in multiple trials are planned, especially in relation to clinicopathological parameters.
Multiregion sampling suggests that colorectal cancer can be highly heterogeneous between tumor regions, with prognostic consequences [4]. In that study, 40% of tumors were found to have CMS variability, typically mixtures of CMS1/3, CMS2/4 or CMS1/4, and usually driven by stroma. Although the iCMS was not developed with the intent of being robust to heterogeneity, it demonstrated a much higher stability across multiregional sampling, with a concordance estimated at 82% (Fig. 4b in [4]). The data presented here do not allow us to evaluate the stability of the iCMS across regions or metastases, but the SSC DQ iCMS distance metric may allow some separation of iCMS mixtures across a continuum, as shown in synthetic data (Figure S4).
Interestingly, tumors with transcriptional similarities to both iCMS2 and iCMS3 (“mixed”) seemed to have unfavorable prognosis in our data, an observation that could point to either a complex ecosystem containing both components, or significant plasticity over time [15]. Indeed, some tumors contain both LGR5+ stem-like and LGR5- regenerative components, which may be driven by reversible epigenetic changes. These components may harmoniously and stably coexist in a mutually beneficial manner. Alternatively, cells may fluctuate between these two states, with the existence of “intermediate” cells arguing in favor of this hypothesis [16]. It seems intuitively likely that this plasticity would enable treatment resistance, but further data, ideally before and after various treatments, are necessary. The formal definition of a “mixed” class of tumors could be an opportunity for further research.
It is unclear how discordant or uncertain iCMS calls could influence patient management, for example if different tumor regions are assigned to different iCMS classes with potentially different treatment sensitivities. The decision to apply a confidence cutoff in a clinical scenario will have to depend on the relative cost of misclassification and the cost of a non-informative result in that specific scenario. This can be configured by the user of the SSC by specifying the minimum distance parameter of the iCMS.DQ function. In practice, the classifier was accurate even for nearest calls, and their use can also be considered for most applications. We hope to that more evidence from multi-region and multi-sample cohorts will demonstrate the stability of the iCMS over time and across lesions.
Several studies suggest that bulk CMS predicts response to therapy. A retrospective analysis of the COIN and PICCOLO trials suggests that an oxaliplatin backbone with anti-EGFR therapy is not effective in CMS1 and CMS4 tumors [3]. In data from the PANAMA trial, CMS2/4 benefited from panitumumab and FU/FA maintenance, in contrast with CMS1/3 tumors [17]. An analysis of three clinical cohorts showed higher response rates to cetuximab combination therapies in CMS2 tumors, even including right-sided CMS2 tumors [18]. Taken together, these results suggest that bulk gene expression profiling may be used to predict treatment response. Compared with bulk CMS, the use of single-cell data should better capture the state of the important epithelial component in iCMS. Indeed, previously published data from cell lines suggest that the iCMS may predict differential sensitivity to commonly used regimens, such as FOLFOX and FOLFIRI, in vitro. Even though clinical data are currently lacking, we hope that the availability of a robust single-sample classifier will facilitate adoption and further exploration of this important question.
The iCMS SSC has limitations. Our validation cohorts were retrospective and may not capture the full spectrum of pre-analytical and analytical methods that could arise in clinical use. The iCMS classification is explicitly focused on the epithelial component and does not account for stroma or immune infiltration. Despite a high proportion of confident calls in the test data, the classifier may be unable to assign a confident iCMS class, possibly because of intratumor heterogeneity (“mixed” tumors) or technical variation. Formal treatment interaction tests in data from interventional randomized clinical trials are necessary to evaluate the ability of the classifier to predict response to therapy.

Conclusions
Further research in colorectal cancer will be facilitated by the existence of a standardized software package for the estimation of epithelial intrinsic molecular subtypes. The iCMS was derived in a non-supervised manner and its prognostic and predictive associations are not fully explored. We hope that by providing a package that is specifically tailored for single-sample analyses, more applications will be made possible. Important open questions remain, notably the ability of the iCMS SSC to predict response to therapy and its performance (and significance) in highly heterogeneous tumors, a topic that we are planning to explore in future studies. We hope that the SSC will be applied retrospectively to trial data, especially in early trials with new compounds, but also in randomized trials to establish true predictive value. Ultimately, the SSC could be used for trial stratification, and clinical decision support if the usual regulatory requirements (platform, laboratory-developed test validation etc) are met.
Exciting new research should complement the epithelial iCMS axis with information from other tissue components, such as the tumor stroma and immune infiltrates. The combination of epithelial, stroma and immune classifiers should be able to better recapitulate tumor biology and presents exciting avenues for further research.

Electronic supplementary material
Below is the link to the electronic supplementary material.