De novo transcriptomic subtyping of colorectal cancer liver metastases in the context of tumor heterogeneity

Background Gene expression-based subtyping has the potential to form a new paradigm for stratified treatment of colorectal cancer. However, current frameworks are based on the transcriptomic profiles of primary tumors, and metastatic heterogeneity is a challenge. Here we aimed to develop a de novo metastasis-oriented framework. Methods In total, 829 transcriptomic profiles from patients with colorectal cancer were analyzed, including primary tumors, liver metastases, and non-malignant liver samples. High-resolution microarray gene expression profiling was performed of 283 liver metastases from 171 patients treated by hepatic resection, including multiregional and/or multi-metastatic samples from each of 47 patients. A single randomly selected liver metastasis sample from each patient was used for unsupervised subtype discovery by nonnegative matrix factorization, and a random forest prediction model was trained to classify multi-metastatic samples, as well as liver metastases from two independent series of 308 additional patients. Results Initial comparisons with non-malignant liver samples and primary colorectal tumors showed a highly variable degree of influence from the liver microenvironment in metastases, which contributed to inter-metastatic transcriptomic heterogeneity, but did not define subtype distinctions. The de novo liver metastasis subtype (LMS) framework recapitulated the main distinction between epithelial-like and mesenchymal-like tumors, with a strong immune and stromal component only in the latter. We also identified biologically distinct epithelial-like subtypes originating from different progenitor cell types. LMS1 metastases had several transcriptomic features of cancer aggressiveness, including secretory progenitor cell origin, oncogenic addictions, and microsatellite instability in a microsatellite stable background, as well as frequent RAS/TP53 co-mutations. The poor-prognostic association of LMS1 metastases was independent of mutation status, clinicopathological variables, and current subtyping frameworks (consensus molecular subtypes and colorectal cancer intrinsic subtypes). LMS1 was also the least heterogeneous subtype in comparisons of multiple metastases per patient, and tumor heterogeneity did not confound the prognostic value of LMS1. Conclusions We report the first large study of multi-metastatic gene expression profiling of colorectal cancer. The new metastasis-oriented subtyping framework showed potential for clinically relevant transcriptomic classification in the context of metastatic heterogeneity, and an LMS1 mini-classifier was constructed to facilitate prognostic stratification and further clinical testing. Supplementary Information The online version contains supplementary material available at 10.1186/s13073-021-00956-1.


Fig. S1. Overview of study material and analyses.
The diagram indicates sample inclusion in the main analyses presented in the study. The study workflow is presented from top to bottom.

Fig. S2 The impact of batch correction on sample type comparisons.
In-house gene expression data were generated in two different batches (batch 1: primary CRCs and cell lines; batch 2: CRLMs, normal liver samples, and PDOs). The plots from PCA shown here are for batch corrected data, and correspond to the plots generated without batch correction in main Fig. 1a-b. PCA was performed on the 1000 genes with highest variation across the batch corrected dataset, and the liver score was calculated by GSVA of a set of genes highly expressed in the liver.  CRLMs were classified according to CMS using our classifier adapted to the liver setting (CMScaller v2.0). Among single CRLMs (randomly selected) from each of the 169 patients, 129 were confidently classified. The pie chart shows the proportion of samples in each CMS group, and the bar plots indicate the proportion of each subtype according to exposure to neoadjuvant chemotherapy (treatment status yes versus no).

Fig. S5. GSEA of the epithelial versus mesenchymal subtype from K=2 factorization.
The plot shows the top 20 enriched pathways in the two subtypes as indicated by colors and ranked according to p-values from GSEA.

a-b
The same plots as Fig. 1a and 1b, respectively, with CRLM samples colored according to LMS. c Distribution of the "liver scores" among the LMS groups indicated no influence of the proportion of hepatocyte signals on de novo subtyping. Wilcoxon test p-values were calculated using LMS1 as reference group, analyzing single samples from each patient (n=169).

Fig. S7. Selected single-sample GSVA scores across the LMS groups.
a MSI-like signature score for each sample (one randomly selected per patient) across the LMS groups. Red asterisk denotes the single sample with confirmed MSI+ status. b Cytotoxic T cell and MSI-like scores plotted by MSI-status in primary CRCs illustrate the relationship between MSI and cytotoxic T cell infiltration in the primary setting. In contrast, the correlation between the MSI-like score and cytotoxic T cell score in CRLMs (mostly MSS) is weak. The linear smooth is represented by the red dashed line. The ρ value represents Spearman's correlation. c KRAS-addiction signature scores in CRLMs with confirmed KRAS mutation plotted according to LMS. Wilcoxon test p-values were calculated with LMS1 as reference group.

Fig. S8. GSEA in CRLMs with RAS/TP53 co-mutations.
Oncogenic pathway enrichments were similar to results from GSEA of all samples. The color intensities represent p-values from comparison of one subtype against all others.

Fig. S9. Kaplan-Meier plots of 5-year CSS according to LMS and translated CMS subtypes.
Kaplan-Meier plots of 5-year CSS stratified by a individual LMS subtype, b LMS1 versus LMS2-5 grouped, and c according to both LMS1 and translated CMS subtypes. Patients with R2 resections in the liver were excluded. Log-rank test p-values are denoted (FDR corrected by the Benjamini-Hochberg in a).

Fig. S10. Kaplan-Meier plots of 5-year OS and CSS according to epithelial and mesenchymal subtypes.
No difference in patient survival was observed when tumors were classified as epithelial or mesenchymal (result from NMF classification at K=2). Long-rank test p-values are denoted. Numbers at risk in each class are shown in the plot with the corresponding color. Patients with R2 resections in the liver were excluded.

Fig. S11. Kaplan-Meier plots of 5-year OS and CSS according to LMS and TP53/RAS co-mutations.
a LMS1 was associated with a poor patient outcome compared to LMS2-5 when analyzing only patients with R0/R1 resections (excluding both patients with R2 resection in the liver, and patients with extra hepatic disease, totally 42 patients). Long-rank test p-values are denoted. b Patient stratification according to both LMS1 versus LMS2-5 and TP53/RAS co-mutation versus no co-mutation showed that there was no significant difference for LMS1 tumors with and without co-mutations. Furthermore, LMS1 tumors with co-mutations were associated with a significantly worse 5-year OS than LMS2-5, both with and without co-mutations. A significant difference in the 5-year CSS rate was only observed in the comparison of co-mutated LMS1 versus LMS2-5 without co-mutation. Numbers at risk in each class are shown in the plot, with corresponding colors. Pairwise FDR-adjusted P-values from log-rank tests are denoted.  Enrichment patterns from GSEA were concordant between each of the two independent datasets and the in-house material (plot corresponding to Fig. 3a). The color intensities represent p-values from comparison of one subtype against all others. a Subtype distributions from de novo subtyping of the in-house (corresponding to the pie chart shown in main Fig. 2d) and two external CRLM datasets. LMS groups of the two external datasets were assigned based on the NMF cluster with the largest number of sample overlaps with supervised LMS predictions. b Comparisons of sample-wise posterior probabilities from the LMS random forest prediction model (vertical axes) and the de novo NMF subtypes (horizontal axes) in each external dataset. Samples with concordant classification between methods (consensus samples) are colored according to the LMS color scheme, or black, as indicated. Samples that are discordantly classified (subtype called by only one of the methods) are colored gray. Notably, each sample was assigned to the subtype with the maximum score for each of the supervised LMS and unsupervised classifications (no fixed threshold for classification).

Fig. S15. Kaplan-Meier curves of OS and CSS according to LMS1 and tumor heterogeneity.
Only patients with R0/R1 resections in the liver were included for analyses. There was no significant difference in the 5-year survival rates a between patients with homogenous versus heterogeneous LMS classifications in inter-metastatic comparisons; or b between patients homogenously classified with LMS1 in all samples versus patients with heterogeneous classifications including at least one LMS1 lesion/sample. LMS1 had poor prognostic associations independent of tumor heterogeneity, as shown by stratification of all patients according to c homogenous LMS1classification (LMS1 in all samples analyzed) versus LMS2-5 plus heterogeneous LMS1 classification; and d LMS1 in at least one lesion versus LMS2-5 in all lesions. Long-rank test p-values are shown. Numbers at risk in each class are shown in the plot, with the corresponding color. Kaplan-Meier plots of 5-year OS and CSS for patient stratification based on the LMS1 mini-classifier. Patients with R2 resections in the liver are excluded. Numbers at risk in each class are shown in the plot, with corresponding colors.