A consensus molecular subtypes classification strategy for clinical colorectal cancer tissues

Comparing FFPE-RNA isolation and library preparation techniques

To establish an effective classification procedure for FFPE CRC tissue, we optimized RNA isolation and library preparation methods using samples from the discovery and FFPE-RNA application cohorts (Table 1, Fig S1). To this end, we first compared the RNeasy FFPE kit, the ReliaPrep FFPE total RNA miniprep kit, and the AllPrep DNA/RNA FFPE isolation kit. RNA-quality, as determined with the mean percentage of transcripts over 200 bp in length (DV₂₀₀), was worst for the former (DV₂₀₀ = 21.7%), and comparable for the latter two procedures (DV₂₀₀ = 37.8%, DV₂₀₀ = 34.0%, respectively). For further isolations, the AllPrep DNA/RNA FFPE isolation kit was chosen because of efficient co-isolation of DNA and RNA. Next, we determined the most effective library preparation strategy for FFPE RNA-sequencing (FFPE-RNA-seq). Six FFPE CRC samples were sequenced with whole exome sequencing (WES), NEB, FASTSELECT, and TAKARA (see the Materials and Methods section for details). To assess transcript preservation, the transcript integrity number (TIN) scores per sample were calculated per protocol, for the genes from the RF CMSClassifier (Fig S2A) and the Single-Sample Predictor (SSP) (Fig S2B). The TAKARA protocol yielded the highest TIN-scores, especially in highly degraded samples. Also, the largest proportion of genes from both the RF and SSP classifiers could be recovered with the TAKARA protocol. Hence, QIAGEN’s AllPrep DNA/RNA FFPE isolation kit and the TAKARA library preparation protocol were used to sequence additional FFPE samples.

A CMS classification model for FFPE archival tissue

Next, we compared matched microarray-derived FF and RNA-seq FFPE gene-expression profiles from the low-quality and moderate-quality discovery sets (Tables 1, S1, and S2). We found a high correlation in both the moderate-quality (Fig 1A) and low-quality set (Fig S3A). However, gene-expression levels of RNA-seq FFPE samples were generally lower than of microarray FF samples. As immune, epithelium, and stroma marker genes are all important to assign the CMSs, we evaluated TIN-scores for these gene categories. Notably, we found that immune marker genes were more likely to be degraded than epithelium and stroma marker genes (Figs 1B and S3B). The degradation of immune genes in the low-quality, moderate-quality, and FFPE-RNA application cohort appeared to be cell type specific, rather than a universal feature of immune cells (Fig S4A–C).

Second, we compared the TIN-score of genes in pathways associated with the CMS subtypes in the discovery sets (4) and observed differences in RNA-degradation level per pathway in FFPE samples, with lowly and highly degraded pathways relatively equally distributed across CMSs (Figs 1C and S3C). This suggests that the level of transcript degradation is also pathway- and sample-specific, rather than CMS-specific, necessitating the use of well-preserved genes for CMS classification of FFPE samples.

Next, we assessed the degradation level of CMS feature genes from the original CMS classifiers in the moderate-quality discovery set. One-hundred forty-four genes met the criterion of a median TIN-score >20 across samples and were selected for classifier building. The expression of these genes in FF and FFPE samples is strongly correlated (PCC = 0.58, P < 0.001) (Fig 1D). As expected, the TIN-scores for these genes were significantly higher than for all other genes (Fig 1E).

With the selected 144 feature genes a new RF classifier was trained on the training cohort (Table 1), termed the CMSFFPE classifier.

Validation of the CMSFFPE classification model

Next, we applied the CMSFFPE classifier to the FFPE samples of the discovery cohort (Fig S1) and compared the results with those predicted by the RF and SSP CMSClassifiers, with core CMS calls from network clustering of matched FF microarray samples as benchmark (4). The overall accuracy of the CMSFFPE classifier was 63.6% in the low-quality set and 83.3% in the moderate-quality set, while for the original RF CMSClassifier the accuracy was 41.0% and 66.7%, respectively, and for the SSP classifier 54.6% and 66.7%, respectively (Figs 2A and S5A–C). For the CMSFFPE model, the predicted CMS calls from FFPE samples and FF samples in the low-quality and moderate-quality sets were highly concordant (P = 0.026, P = 0.0040, respectively, Binomial Exact test, Fig S5D), with the highest specificity for CMS2 (both 100%), superior to the RF CMSClassifier and the SSP (Fig S5E and F). Importantly, for the RF CMSClassifier, the pairwise concordance was not statistically significant (P = 0.58, P = 0.073, respectively, Binomial Exact test, Fig S5E), underlining the need for an FFPE-curated classifier.

To explore the performance of the new CMSFFPE classifier in a larger and high-quality FF-sample dataset, we used microarray-derived RNA data of FF samples from the Cartes d’Identité des Tumeurs (CIT) cohort (FF-RNA validation cohort, Table 1, Fig S1) (11). As a benchmark, we used core CMS calls from network clustering (4). Subsequently, the CMSFFPE and RF CMSClassifier were applied to predict the CMS subtypes. For both classifiers, the CMS calls were highly concordant with the benchmark calls (Fig 2B). The overall accuracy of the CMSFFPE classifier (89.0%) was similar to the RF CMSClassifier (88.4%) and the SSP classifier (91.1%). Notably, discordant CMS calls between the CMSClassifier and the CMSFFPE classifier manifested as a shift from, respectively, CMS1/3/4 toward CMS2, reflecting the higher accuracy of the CMSFFPE classifier to detect CMS2 (81.5% versus 93.1%). Samples that were discordantly classified showed a significantly lower CMSFFPE probability score, suggesting that these samples had a more CMS-mixed phenotype (Fig S6A–D). Importantly, CMS-dependent differences in relapse-free survival (RFS) could be recapitulated with the new CMSFFPE classifier and were more significant than those identified by the original RF CMSClassifier (Fig 2C). With both the new and the original classifier, CMS1 had the best and CMS4 the worst RFS (Hazard ratio [HR] 2.61, [95% confidence interval (CI) 1.46–4.67]; HR 2.15 [1.29–3.59], respectively).

Application of the CMSFFPE classifier to the FFPE-RNA application cohort

Because previous studies have provided convincing support for differential anti-EGFR sensitivity across the CMSs (17), we evaluated the performance and clinical relevance of the CMSFFPE classifier in an independent anti-EGFR treated FFPE cohort with a matched control group, consisting of 104 RAS/BRAF wt mCRC patients (Table 1, Fig S1). To reflect clinical reality, we obtained either a biopsy or resection specimen of the primary tumor (Table 2). Sixty-two patients were treated with anti-EGFR and these had more distant metastases, which were regularly metachronous, were treated with more lines of systemic treatment and had a worse overall survival (OS) compared with the control group.

As a data-quality check, we assessed the transcriptome mapping and coverage for all samples (Table S3). Most of the samples (90.4%) showed low mapping rates (<50%) (Fig S7A) (25), indicating moderate data quality. Moreover, the fraction of intronic reads was relatively high, ranging from 44.2%–64.5% (Fig S7B). As a consequence of the used library preparation protocol, we found that the read fragments contained large proportions of adapters in the raw sequencing data (Fig S7C), with a median insert size of 72 (Fig S7D). Consistent with the discovery cohort, immune marker genes had the lowest TIN (Fig S7E), and the degradation levels of genes per pathway closely resembled those described for the discovery cohort, despite the different RNA-seq protocols (Fig S7F). In addition, we compared data quality of biopsy versus resections specimens and found a higher quality of biopsy specimens, likely related to a shorter fixation time of biopsy samples (Fig S8A–C). Upon CMSFFPE classification, we identified 4 (3.8%) CMS1 samples, 58 (55.8%) CMS2 samples, 10 (9.6%) CMS3 samples, and 32 (30.8%) CMS4 samples conform the expected distribution in RAS/BRAF wt mCRC patients (26).

Molecular characteristics of the CMSs in the FFPE-RNA application cohort

To confirm that the CMSs in the FFPE-RNA application cohort were molecularly similar to the original CMSs, we performed molecular characterization per CMS (Fig 3). First, differential expression analyses of classifier genes showed proper segregation of the CMSs, confirming robust performance of the CMSFFPE classifier (Fig 3A). Next, we assessed gene set enrichment across the CMSs and could recover all previously reported iCMS features (Fig 3B). These included enrichment for immune activation in CMS1, SRC signaling, cell cycle and DNA repair in CMS2, metabolic processes in CMS3 and mesenchymal differentiation, EMT, and TGF-β signaling in CMS4 (27). Furthermore, the immune cell signature was enriched in CMS1, differentiated enterocytes in CMS2, goblet cells in CMS3, and mesenchymal cells in CMS4.

In addition, we performed mutational characterization with the TruSight Oncology 500 (TSO500) assay and found that patients were most frequently TP53- (75%) and APC-mutated (57%), probably related to the RAS/BRAF wt status and predominant left-sidedness of the cohort (Fig 3C). Two patients with MSI tumors were classified as CMS1 and CMS4. Compared with the other CMSs, CMS2 tumors were enriched for canonical mutations in TP53 (89.1%, P = 0.002), and were relatively more frequently mutated in APC (63.6%, P = 0.312), TP53 and APC combined (58.2%; P = 0.321), and EGFR (21.8%; P = 0.085) (Table S4). Altogether, these findings corroborate robust CMS classification with the CMSFFPE classifier.

Clinical characteristics per CMS in the FFPE-RNA application cohort

Furthermore, we assessed whether clinical associations of the CMSs could be recovered in our FFPE-RNA application cohort. In line with the RAS/BRAF wt status of the cohort, most tumors were left-sided, except for CMS1 (75% right-sided) (Table 3). CMS2 patients nearly all had liver metastases (CMS2 versus rest: P = 0.004; Chi-square test) and a relatively long OS (median OS 42.9 mo [38.5–47.4]), whereas CMS3 tumors showed the worst OS (median OS 20.7 mo [10.7–30.6]) (data not shown). This might be related to the rather large proportion of advanced CMS3 patients with synchronous metastases compared with the other subtypes. Strikingly, CMS4 tumors performed best with an OS comparable with CMS2 cancers (median OS 47.6 mo [37.9–57.3]). We also assessed OS for the anti-EGFR treated and untreated subgroups separately and found a similar pattern in both groups, with CMS2 and CMS4 having a comparable OS (HR 0.86 [0.45–1.66], HR 1.07 [0.39–2.90], respectively) (Fig S9A and B). In univariable and multivariable Cox regression models, the CMSs did not independently predict OS (P = 0.183), although OS patterns resembled known prognostic associations of the subtypes in mCRC (17), with CMS2 having a better OS than the other subtypes (Table S5, Fig S10).

Molecular predictors of anti-EGFR response

As a proof of principle for the feasibility of FFPE-RNA and -DNA based precision medicine, we compared mutations between patients that responded to anti-EGFR therapy and patients that did not respond (Fig 4A). We found a correlation between anti-EGFR response and TP53 mutations (96.2% in responders versus 63.6% in non-responders, P = 0.017) (Table S6). Strikingly, the combination of a TP53 and APC mutation was even stronger associated with anti-EGFR response (P = 0.004).

Because previous literature suggested a benefit of anti-EGFR therapy for CMS2 tumors, we assessed whether the same patterns could be observed in our FFPE-RNA application cohort. Because CMS1 and CMS3 tumors were scarce in the treated group, we compared anti-EGFR response rates of epithelial (CMS2&3) versus immune and mesenchymal (CMS1&4) tumors, and observed a (non-significant) pattern of a better objective response rate (ORR) for epithelial tumors (ORR 62.5% versus 38.9%, respectively) (Fig 4B). Interestingly, two right-sided tumors responded to anti-EGFR therapy, which were both epithelial tumors with co-occurring APC and TP53 mutations. Among non-responders were six right-sided tumors, of which two were epithelial and two had co-occurring APC and TP53 mutations. After univariable and multivariable logistic regression analyses, independent predictors for anti-EGFR therapy responsiveness (ORR) were CMS2&3 tumors (P = 0.038) and co-occurring TP53 and APC mutations (P = 0.007) (Table S7, Fig 4C). This suggests that FFPE-based CMS classification and gene-panel DNA sequencing are robust methods to detect distinct patient groups sensitive to systemic therapies.

Patient series

Five patient cohorts were used for, respectively, selection of well-preserved genes in FFPE samples (discovery cohort – moderate-quality set), classifier training (training cohort), classifier testing (discovery cohort – low- and moderate-quality sets), classifier application on FF-derived RNA data (FF-RNA application cohort), and classifier application on FFPE-derived RNA data (FFPE-RNA application cohort) (Table 1). The datasets were used in multiple stages of the study (Fig S1).

The discovery sets consisted of patients with localized CRC that underwent curative surgery (6) (Table 1). Clinical information, including survival data, was available. For all discovery patients, both FFPE-derived and FF-derived RNA were available, originating from RNA-seq and microarray analysis, respectively. FF-derived core CMS calls from network clustering were used as benchmark (4). Well-preserved genes in FFPE samples were selected, which were used to construct the CMSFFPE classifier.

The cohorts of colon adenocarcinomas (COAD) and rectal adenocarcinomas (READ) from the Cancer Genome Atlas (TCGA) were used for classifier training (training cohort) (2) (Table 1). CMS calls from the original RF CMSClassifier served as benchmark.

As a microarray-derived FF-RNA validation cohort, the CIT cohort was used (Table 1). This cohort consisted of 519 CRC patients that underwent surgery of the primary tumor, with documented RFS and core CMS calls from network clustering (4). The core CMS calls were used as benchmark for comparing the accuracy of the CMSFFPE classifier with the original RF CMSClassifier.

As an FFPE-RNA application cohort, a total of 104 FFPE tumor biopsy or resection samples from unique RAS/BRAF wt mCRC patients diagnosed from January 2015 through December 2019 were requested from the Prospective Dutch CRC cohort study (PLCRC) (46). The sample size gave >95% certainty to detect at least one molecular variant with a prevalence of 3% in the total cohort (power calculation: 0.97ⁿ ≤ 0.05; n = sample size). All patients provided written informed consent for usage of their clinical data and tumor specimens. PLCRC was registered at ClinicalTrials.gov (NCT02070146) and approved by the Medical Research Ethics Committee of the University Medical Centre Utrecht (NL47888.041.14). PLCRC is conducted in accordance with the declaration of Helsinki (46).

Patients were selected based on previous treatment with anti-EGFR (the treatment population), and no previous treatment with anti-EGFR (the control population), and matched on age, sex, ASA-classification, stage, and sidedness. Anonymized clinical data were obtained via the Netherlands Cancer Registry (47). Clinical data consisted of patient characteristics (age, sex, survival), tumor characteristics (stage at diagnosis, metastases, lymph node involvement), and treatment characteristics (surgery, chemotherapy, targeted therapy, documented best response to therapy). Clinical data and survival were last updated in March 2022.

FFPE-RNA isolation

We compared three different FFPE-RNA isolation procedures on RNA quality, as estimated with the DV₂₀₀, determined by the 4200 TapeStation System (Agilent). We compared the RNeasy FFPE kit (QIAGEN), the ReliaPrep FFPE total RNA miniprep kit (Promega) and the AllPrep DNA/RNA FFPE isolation kit (QIAGEN). For the FFPE-RNA application cohort, total RNA and genomic DNA were isolated from FFPE samples with QIAGEN’s QIAcube classic (cat. no. 9001293; QIAGEN). A DV₂₀₀>15% was deemed sufficient for RNA-seq.

The transcript integrity number

We assessed the transcript preservation levels of our FFPE datasets by calculating the TIN (28). The TIN reflects the percentage of a transcript with uniform read coverage. The TIN-score was calculated for each sample using RSeQC’s tin.py (v.5.0.1) for the annotated transcripts (48).

RNA-sequencing and targeted exome sequencing

To account for variations in FFPE-RNA-sequencing methods, two different procedures were used for FFPE samples of the discovery cohort, dividing the cohort into a low-quality and moderate-quality set, as objectified with the median TIN-score (Tables 1, S1, and S2).

Before sequencing the FFPE-RNA application cohort, four library preparation strategies were compared in six CRC FFPE samples, to determine the impact of library preparation strategy on data quality. We compared (A) whole exome sequencing, (B) NEB Ultra II Directional library preparation with rRNA depletion (NEB), (C) NEB cDNA module with rRNA depletion, Agilent SureSelect All Exon V7 library preparation (FASTSELECT) and (D) SMARTer Stranded Total RNA-seq Kit v3 – Pico from TakaraBio (TAKARA). The library preparation technique yielding the optimal RNA-data quality was applied to the full cohort.

The TSO500 High-Throughput assay was performed on genomic DNA of the FFPE-RNA application cohort as instructed by IlluminaPaired-end next-generation sequencing was performed on the NovaSeq6000 by GenomeScan. The TSO500 detects somatic variants across 523 cancer-relevant genes, as well as tumor mutational burden (TMB) and mismatch repair (MMR) status. Classification of TMB-high or TMB-low was based on a cut-off value of 10 mutations per megabase. MMR status was determined as a percentage based on the analyses of 130 MSI marker sites. A cut-off of 15% was chosen to distinguish MSI-high from MSI-low.

RNA-data processing and gene-expression quantification

For the discovery cohort, fragments were mapped to the human genome (GRCh38) using STAR v2.7.8a with GENCODE v33 and default parameters (49). Only uniquely mapped fragments were considered and transcripts per million (TPM) were calculated using RSEM v1.3.3 (50). This methodology was similar to TPM quantification performed on TCGA data, which was used for classifier training. Hence, we excluded bias because of differences in TPM quantification methods in the datasets used for feature selection and classifier training.

For the FFPE-RNA application cohort, the RNA-seq data were analyzed using the Cogent next-generation sequencing Analysis Pipeline (CogentAP 2.0), as suggested by TakaraBio (https://www.takarabio.com/products/next-generation-sequencing/bioinformatics-tools/cogent-ngs-analysis-pipeline). In brief, we removed sequencing adapters and trimmed the RNA-seq reads using cutadapt (51) (v3.4, settings: -m 15 --trim-n --max-n 0.7 -q 20). Next, the retained RNA-seq reads were mapped to the human reference genome (GRCh38) and the reads were annotated with GENCODE v33 using STAR (49). The number of uniquely mapped reads per gene per sample was calculated using featureCounts (v2.0.0) (52). Of note, RSEM and featureCounts to quantify TPM were previously shown to be highly concordant (53). Gene count matrices were imported into R (v4.0.3) and converted into TPM with the genomeutils package (https://github.com/ssarda/genomeutils). The Picard software (https://broadinstitute.github.io/picard/) (v2.27.3) was used to estimate the insert sizes of reads.

Redesigning the CMS classifier

The CMSFFPE classifier was developed with the RandomForest R package (54). The gene selection and model building steps are as follows:

• (i) Gene selection: signature genes of the CMSClassifier (4) and CMScaller (55) with a median TIN-score >20 (56) in the moderate-quality set were kept as features for the CMSFFPE classifier.

• (ii) Training set: the classifier was trained on the TCGA-CRC cohort with core CMS calls from network clustering (4) as benchmark. TPM values were downloaded with the R package “TCGAbiolinks” (57) in “harmonized” mode.

• (iii) Model building: to account for variations in gene expression caused by distinct experimental protocols, we first performed z-normalization on the gene-expression profiles, and then derived binary gene pairs from the standardized expression values. This methodology has been previously used for ovarian cancer subtyping (58).

• (iv) Validation in FFPE and FF-RNA datasets: to evaluate the performance of the CMSFFPE classifier, we assessed the concordance between the CMS calls predicted from FFPE samples and paired FF samples in the discovery cohort, and compared them with the calls derived from the original CMSClassifier. Next, we applied the CMSFFPE classifier and the CMSClassifier to the FF-sample data of the CIT cohort (GSE39582) (11).

• (v) Application in an independent FFPE-RNA dataset: lastly, we applied the CMSFFPE classifier to the FFPE-RNA application cohort. Molecular features per CMS were assessed by analyzing differential expression of genes (DEG) per CMS compared with the others, and visualized with a heatmap, using the R packages “edgeR” and “pheatmap.” Gene set enrichment analyses for samples with read counts >10⁷ were performed and visualized using the R packages “fgsea” and “pheatmap.”

Analyses of TSO500 mutation data

Fastq data were summarized in standardized combined variant output reports, using the local app v2.2.0 docker container (Illumina) executed on a local HP Unix server via Singularity. Combined variant output reports were parsed to select variants with an amino-acid change, a frameshift, splice variants and start- and stop-events. To account for differences in tumor sample purity, tumor cell percentages, as assessed by a certified pathologist, were used to correct variant allele frequencies. To select for driver variants, the following filter criteria were used: corrected variant allele frequency >0.10 for oncogenes and >0.40 for tumor suppressor genes; presence of a variant in the population ≤10⁻⁵. To avoid inclusion of assay-specific sequence artifacts, an in-house reference database was set up that included more than 450 TSO500 profiled tumor samples. This reference was used to filter out those variants that occurred in more than 50 percent of the tumors. To avoid missing important somatic mutations, variants that were seen ≥1 cases of the TCGA (from cBioportal) (COAD and READ), were maintained. Also, somatic variants reported in ≥10 cases of the Memorial Sloan Kettering Cancer Center (from cBioportal) (pan cancer) were maintained.

Cohort and sample overviews and visualizations were implemented in the R2 Genomics Analysis and Visualization Platform (https://r2.amc.nl), together with sample annotations. Furthermore, oncoplots were generated to visualize mutations per CMS and per anti-EGFR response group, using the R package “ComplexHeatmaps” from BioConductor (59).

Associations of subtypes with clinical data

Analyses of clinical data were performed separately from analyses of molecular data, by different researchers. SPSS statistical software, version 28.0.1.1, and R version 4.0.5, were used to perform analyses into clinical correlates of molecular characteristics per tumor. Non-parametric tests were carried out to analyze correlations between continuous clinical data and mutations or CMSs. The Chi-square test, the Fisher’s exact test or the Fisher-Freeman-Halton Exact test was used to associate categorical data with molecular characteristics of tumors. Univariable and multivariable logistic regression analyses were performed to identify most significant determinants of anti-EGFR response in the treatment group. Response was defined as the ORR: partial or complete response as registered by the treating physician.

Survival analyses were carried out using the survival package in R. OS was defined as time from diagnosis until death of any cause. Survival curves for molecular and clinical subgroups were generated with the Kaplan-Meier method and compared with the log-rank test. Cox proportional hazards modeling for OS was performed, after confirming the assumption of proportional hazards. Univariable and multivariable analyses were carried out to identify prognostic determinants for OS.

Most significant and most clinically relevant variables from univariable analyses were included in multivariable analyses. Variables with missing values were excluded from multivariable analyses. P-values were two-sided with an alpha < 0.05. Multiple testing correction was performed for analyses with sequencing data.

Author Contributions

• TR de Back: conceptualization, data curation, software, formal analysis, validation, investigation, visualization, methodology, and writing—original draft.

• T Wu: conceptualization, data curation, software, formal analysis, validation, investigation, visualization, methodology, and writing—original draft.

• PJM Schafrat: conceptualization and writing—review and editing.

• S ten Hoorn: conceptualization, methodology, and writing—review and editing.

• M Tan: formal analysis and investigation.

• L He: formal analysis and investigation.

• SR van Hooff: data curation, formal analysis, investigation, and visualization.

• J Koster: data curation, formal analysis, investigation, and visualization.

• LE Nijman: investigation.

• GR Vink: conceptualization and writing—review and editing.

• IJ Beumer: conceptualization and writing—review and editing.

• CC Elbers: conceptualization, supervision, funding acquisition, project administration, and writing—review and editing.

• KJ Lenos: conceptualization, supervision, funding acquisition, and writing—review and editing.

• DW Sommeijer: conceptualization, resources, supervision, funding acquisition, validation, methodology, project administration, and writing—review and editing.

• X Wang: conceptualization, resources, supervision, funding acquisition, validation, methodology, project administration, and writing—review and editing.

• L Vermeulen: conceptualization, resources, supervision, funding acquisition, validation, methodology, project administration, and writing—review and editing.

Conflict of Interest Statement

L Vermeulen received consultancy fees from Bayer, MSD, Genentech, Servier, Roche, Novartis, and Pierre Fabre, but these had no relation to the content of this publication. L Vermeulen is currently an employee of Genentech Inc. and shareholder of Roche. GR Vink received consultancy fees from BMS, Merck, Servier, Personal Genome Diagnostics, Bayer, Sirtex, Pierre Fabre, Lilly, and Delfi Diagnostics, but these had no relation to the content of this publication. The other authors declare no potential conflicts of interest.