Multiregional transcriptomics identifies congruent consensus subtypes with prognostic value beyond tumor heterogeneity of colorectal cancer

Abstract

Intra-tumor heterogeneity compromises the clinical value of transcriptomic classifications of colorectal cancer. We investigated the prognostic effect of transcriptomic heterogeneity and the potential for classifications less vulnerable to heterogeneity in a single-hospital series of 1093 tumor samples from 692 patients, including multiregional samples from 98 primary tumors and 35 primary-metastasis sets. We show that intra-tumor heterogeneity of the consensus molecular subtypes (CMS) is frequent and has poor-prognostic associations independently of tumor microenvironment markers. Multiregional transcriptomics uncover cancer cell-intrinsic and low-heterogeneity signals that recapitulate the intrinsic CMSs proposed by single-cell sequencing. Further subclassification identifies congruent CMSs that explain a larger proportion of variation in patient survival than intra-tumor heterogeneity. Plasticity is indicated by discordant intrinsic phenotypes of matched primary and metastatic tumors. We conclude that multiregional sampling reconciles the prognostic power of tumor classifications from single-cell and bulk transcriptomics in the context of intra-tumor heterogeneity, and phenotypic plasticity challenges the reconciliation of primary and metastatic subtypes.

Fig. 1. Landscape of CMS heterogeneity among multiregional samples of primary CRCs.

a CMS classification of 2–4 multiregional samples from each of 98 primary tumors ordered according to intra-tumor classification heterogeneity. Each column represents one tumor. Spatially separated samples are annotated t1–t4 and colored according to the consensus molecular subtypes (CMS). MSI, RAS/BRAF and Location indicate the microsatellite instability status, mutation status (for BRAF^V600E, KRAS and NRAS) and location of each tumor in the large bowel. b Box plot of the general transcriptomic heterogeneity of each tumor, estimated as the maximum Euclidean distance of PC1–PC3 between any pairs of samples per tumor and plotted according to CMS heterogeneity. The center line of boxes represents the median, boxes represent the interquartile range, and whiskers represent 1.5× the interquartile range above the 75th percentile (maxima) or below the 25th percentile (minima). p value is two-sided and from Welch’s t-test. Source data are provided as a Source Data file. c Frequency of intra-tumor heterogeneity of each CMS class. The odds ratio and 95% confidence interval for enrichment of heterogeneous tumors are showed above each class. Source data are provided as a Source Data file. d Bar plot of gene sets sorted according to p values from gene set enrichment analysis of tumors with heterogeneous versus homogeneous CMS classifications (one random sample per tumor). Enrichment in heterogeneous or homogeneous tumors is indicated by the color code. Source data are provided as a Source Data file.

Fig. 2. Prognostic value of intra-tumor CMS heterogeneity.

a Relapse-free survival according to intra-tumor heterogeneity of the consensus molecular subtypes (CMS) among patients treated by complete resection of stage I–III colorectal cancer. Patients with synchronous tumors, pre-surgical radiation treatment and undetermined CMS heterogeneity were excluded from all analyses. b Explained variation in 5-year relapse-free survival by each variable in a multivariable Cox proportional hazards model, estimated as the percentage of the full model. c Survival according to intra-tumor CMS heterogeneity and stratified by CMS4 classification. Hazard ratios (HR) and 95% confidence intervals (CI) are from Cox proportional hazards analyses and p values are two-sided and from Wald tests. CAF cancer-associated fibroblasts, CTL cytotoxic lymphocytes, het heterogeneous, hom homogeneous, TN stage tumor-node stage.

Fig. 3. Gene categories and enrichments according to intra-tumor heterogeneity.

a Upper part shows the proportion of genes (filtered to include only genes with inter-sample 10–90th percentile expression range >1) categorized as either ITH-low (n = 1540), ITH-intermediate (n = 1549) or ITH-high (n = 149) in the multiregional primary tumor set (n = 286 samples). Lower part shows the density distribution of the corresponding ITH-scores, with dashed lines indicating thresholds for the three gene categories. The forest plot shows the ITH-score of genes designated as colorectal cancer-associated in the Cancer Gene Census (n = 13 genes passing the filter for the inter-sample expression range). The genes are sorted according to involvement in pathways over-represented among ITH-low genes (Wikipathways cancer; Supplementary Fig. 14). Source data are provided as a Source Data file. b Density plots (upper: all gene sets, n = 54) and bar plots (lower: 15 top-ranked gene sets) of Spearman’s correlation coefficients (absolute values) between PC1 or PC2 of ITH-high (right) or ITH-low (left) genes and single-sample enrichment scores of significantly correlated gene sets (p < 0.05). Source data are provided as a Source Data file. c Number of gene sets (of totally 54) with significant Spearman’s correlations (p < 0.05) to PC1 and PC2 of ITH-high and ITH-low genes. Source data are provided as a Source Data file. d Distribution of ITH-scores among genes in selected signatures. The signatures were grouped into six categories, and up to three of the top-ranked gene sets according to the correlation analyses in part b are plotted per category, in addition to the WNT/β-catenin signature. Source data are provided as a Source Data file. CIN chromosomal instability, EMT epithelial-mesenchymal transition, ITH intra-tumor heterogeneity, MSI microsatellite instability, MSS microsatellite stable, PC1 and PC2 principal components 1 and 2, ρ Spearman’s correlation coefficient.

Fig. 4. Classification of primary tumors and liver metastases based on ITH-low and cancer cell-intrinsic genes.

a Alluvial diagrams of classification concordances between iCMS and the k2 clusters identified based on ITH-low genes among primary colorectal tumors and liver metastases. The sample overlap between classes in iCMS and k2 is indicated relative to the total number per subtype. Source data are provided as a Source Data file. b Proportion of primary tumors with homogeneous and heterogeneous intra-tumor classifications of multiregional samples (n = 286 samples) according to the indicated frameworks (k2 and k4 are based on the ITH-low genes). Proportion of patients with inter-metastatic heterogeneity among liver lesions (2–7 distinct lesions per patient, total n = 143 metastases) according to iCMS and the k2 clusters is plotted to the right. Source data are provided as a Source Data file. c iCMS classifications of matched primary tumors and liver metastases from 35 patients (n = 179 samples). Each column represents one patient. For primary tumors, each square represents one multiregional sample numbered with lower case t. For liver metastases, each square represents one tumor numbered with upper case T and separated by consecutive resections, and diagonal lines indicate multiregional samples. Indicated for each patient is the platform used for gene expression analysis, diagnosis with synchronous versus metachronous metastases and exposure to chemotherapy prior to sampling. cCMS congruent consensus molecular subtypes, CRIS colorectal cancer intrinsic subtypes, iCMS intrinsic consensus molecular subtypes, ITH intra-tumor heterogeneity, k2 and k4 factorization ranks 2 and 4 from non-negative matrix factorization.

Fig. 5. Prognostic value of the proposed congruent CMS framework.

a Proportion of the four cCMS classes among primary colorectal tumors (n = 704 samples from 516 tumors). Source data are provided as a Source Data file. b Alluvial diagrams of classification concordances between k2 and cCMS (left; representing k2 and k4 classifications based on ITH-low genes), as well as cCMS and CMS (right). The sample overlap between classes is indicated relative to the total number per subtype and illustrated by dashed lines. Source data are provided as a Source Data file. c Proportion of MSI tumors across the cCMS classes. Source data are provided as a Source Data file. d Relapse-free survival according to cCMS among patients treated by complete resection of stage I–III CRC (n = 398). Patients with heterogeneous intra-tumor cCMS classifications, synchronous tumors, or pre-surgical radiation treatment were excluded from analyses. Hazard ratios (HR) and 95% confidence intervals (CI) are from Cox proportional hazards analyses and p values are two-sided and from Wald tests. e Explained variation in 5-year relapse-free survival (n = 398 patients) by each variable in a multivariable Cox proportional hazards model, estimated as the percentage of the full model. CAF cancer-associated fibroblasts, cCMS congruent consensus molecular subtypes, CTL cytotoxic lymphocytes, iCMS intrinsic consensus molecular subtypes, ITH intra-tumor heterogeneity, k2 and k4 factorization ranks 2 and 4 from non-negative matrix factorization, MSI microsatellite instability, MSS microsatellite stable, TN stage tumor-node stage.

Fig. 6. Classification of external primary tumor series based on ITH-low genes.

Alluvial diagrams of classification concordances (a) between iCMS and the ITH-low k2 clusters among tumors in the GSE39582 series, (b) the ITH-low k2 and k5 clusters and the original CMS in GSE39582 (top), as well as the iCMS, ITH-low k6 and original CMS in TCGA (bottom). The sample overlap is indicated relative to the total number per subtype. Source data are provided as a Source Data file. c Pie charts of the proportion of tumors in each of the ITH-low k5 and k6 clusters in GSE39582 and TCGA, respectively, bar charts of the proportion of MSI tumors in each cluster, and (d) heat maps of p values from gene set enrichment analyses (log10-scale; red: positive; blue: negative). Gene sets were from the custom collection and selected to include the same as in the corresponding analysis of the in-house tumor series (Supplementary Fig. 20). Source data are provided as a Source Data file. e Five-year relapse-free survival according to the ITH-low k5 clusters in patients with stage I–III CRC in GSE39582. Hazard ratios (HR) and 95% confidence intervals (CI) are from Cox proportional hazards analyses and p values are two-sided and from Wald tests. The two-sided log-rank p value across subtypes is also given. f Common pathway enrichments in the KEGG database of differentially expressed genes between the two CMS2-corresponding subtypes (NMF2 and NMF2.5) from ITH-low k5 clustering of the in-house and GSE39582 series, and ITH-low k6 clustering of the TCGA series. The eight pathways with highest max enrichment score among the common pathways of the top-50 most significantly enriched in each tumor series were selected for plotting. Dot sizes indicate the significance level (Benjamini–Hochberg adjusted p value on log10-scale). Source data are provided as a Source Data file. cCMS congruent consensus molecular subtypes, EMT epithelial-mesenchymal transition, iCMS intrinsic consensus molecular subtypes, ITH intra-tumor heterogeneity, k2 and k4 factorization ranks 2 and 4 from NMF, MSI microsatellite instability, MSS microsatellite stable, NMF non-negative matrix factorization, TCGA The Cancer Genome Atlas.