Molecular dissection of colorectal cancer in pre-clinical models identifies biomarkers predicting sensitivity to EGFR inhibitors

Abstract

Colorectal carcinoma represents a heterogeneous entity, with only a fraction of the tumours responding to available therapies, requiring a better molecular understanding of the disease in precision oncology. To address this challenge, the OncoTrack consortium recruited 106 CRC patients (stages I-IV) and developed a pre-clinical platform generating a compendium of drug sensitivity data totalling >4,000 assays testing 16 clinical drugs on patient-derived in vivo and in vitro models. This large biobank of 106 tumours, 35 organoids and 59 xenografts, with extensive omics data comparing donor tumours and derived models provides a resource for advancing our understanding of CRC. Models recapitulate many of the genetic and transcriptomic features of the donors, but defined less complex molecular sub-groups because of the loss of human stroma. Linking molecular profiles with drug sensitivity patterns identifies novel biomarkers, including a signature outperforming RAS/RAF mutations in predicting sensitivity to the EGFR inhibitor cetuximab.

Figure 1. Experimental design of the OncoTrack (OT) study.

Resected CRC patient tumours were fragmented and sampled for fuelling the sequencing and the establishment of in vivo PDX and in vitro PDO models. Untreated original tumours, PDX and PDO samples were analysed by WGS, WES and RNAseq for correlating the molecular information with drug sensitivity patterns. In addition, the epigenomes of the original tumours were analysed. The OT pre-clinical platform treated both model systems with therapeutic compounds representing the standard of care and/or addressing major pathways relevant in CRC.

Figure 2. Genomic landscape of the OT patient and model cohorts.

(a) Examples of gene fusions, either deleterious (ME2-SMAD4 and APC-REEP5, blue) or activating (TRIM24-BRAF and FDFT1-FZD3, red); schematics display chromosomal location, fusion partners, exon structure (blue, red and grey filled boxes), fusion breakpoint (red line), validation by Sanger sequencing and RNA read coverage. (b) Landscape of the most recurrent somatic alterations in 101 OT primary and metastasis samples (tumour purity of ≥20%) from 96 patients compared with the mutation pattern in 228 primary CRC tumours from TCGA (ref. 25) (SNV+Indel). Genes are grouped according to biological pathways. MSI and hypermutated samples are shaded in grey. Tumour stages are coded in grey shades on top of the matrix. Metastasis and primary samples are indicated respectively in purple and green boxes on top of the matrix. Alteration types are colour-coded as indicated. Obtained PDX and PDO models are depicted in blue and yellow, respectively at the bottom of the matrix. Darker shades correspond to the models shown in Fig. 2c. (c) Fraction of damaging and expressed somatic SNVs/Indels found discordant between original patient tumours and matched models, comparing 56 patient samples (tumour purity ≥40%) and their corresponding ex vivo models (37 PDX and 37 PDO), as well as five PDX-PDO siblings without matching patient samples. Dot colours indicate the following: patient versus PDX (blue), patient versus PDO (yellow) or PDX versus PDO (red). Samples described in Fig. 3b are marked in bold. MSI and hypermutated samples are marked with an asterisk. Panels from left to right: patients with an established PDX only; patients with an established PDO only; patients with established PDX and PDO; PDX/PDO siblings without a corresponding sequenced patient tumour. In the cases of tumours 278 and 302, the suffix ‘X' represents patient-derived PDX and ‘CX' represents PDO-derived PDX.

Figure 3. Molecular discordances between original tumours and derived models in CRC-relevant genes.

(a) Examples of clonality analysis based on the sciClone algorithm, for samples 327_T (MSI—40% tumour purity) (upper panel) and 323_T (MSS—60% tumour purity) (lower panel). The plots display the mutations in diploid regions clustered by their allele frequencies in: original tumour versus PDX (left), original tumour versus PDO (middle) and PDX versus PDO (right). Individual mutation clusters were shown by different colours where grey indicated the mutations found in common between the compared samples. The PDO derived from 327_T displayed private EGFR, MLL2 and PIK3CA mutations in cluster 2. The PDX derived from 323_T displayed different clones harbouring private mutations in PIK3CA. (b) Somatic clonal mutations discordant between patients/PDX (left), patients/PDO (middle) and patients/PDX/PDO (right). Genes found mutated only in either tumour or model (or in only one model) are marked by red squares. Genes with different mutations in tumours and models are marked by orange squares. Mutations validated by targeted sequencing are indicated with asterisk. MSI and hypermutated samples are shaded in grey: Only cancer-relevant genes are shown, selected from TCGA, PAN or OT recurrently mutated genes resources (see methods). Patients with a purity of ≥40% are shown.

Figure 4. Expression profiling of the OT patient (tumour purity ≥40%) and model cohorts.

(a) Unsupervised hierarchical clustering of the mean patterns of 38 CRC-related signatures in 90 patients, 53 PDX and 33 PDO models. Information on molecular groups is indicated as follows from top to bottom under the dendrogram. Colour labels of the three OT patient tumour groups, NMF classes for the models, and CMS groups (RF: upper row, SPP: bottom row) are indicated in the caption on the right side. MSI- and hypermutated samples are marked with black and blue boxes, respectively. (b) Sankey plots showing the correspondence between the molecular groups in the patient tumours (ECM/EMT, ASCL2/MYC, Entero/Goblets) and the two NMF groups in PDX and PDO models. (c) Sankey plots showing the correspondence between the NMF molecular groups between sibling pairs of PDX and PDO models. (d) Analysis of WNT-signalling reporter assays on 151_MET PDO cultures. Left: cell cultures showing the GFP green fluorescent signal revealing WNT signalling and FACS-sorting of the corresponding GFP-high, -low and -negative cell fractions. Scale bar is 100 μm. Right: Heatmap representing differentially expressed genes (n=404) between the GFP-low, GFP-high, GFP-negative cell fractions and unsorted cells. Genes were filtered by fold change (FC) and difference of RPKM (DR, |DR|≥1, |log₂(FC)|≥2). Key enriched processes and P-values in WNT-negative and -high fractions are indicated on the right side.

Figure 5. Comparative gene expression profiling between PDO and PDX models.

(a) Differentially expressed canonical pathway gene signatures. Positive or negative normalized enrichment scores (NES) calculated for PDX versus PDO for each given pathway are indicated in purple or green shades, respectively. Asterisks indicate significant enrichment with FDR<0.005. The gene set enrichment was determined using a pre-ranked GSEA. Genes were ranked according to their expression fold change and expression difference between the corresponding model pairs. (b) Differential expression of gene signatures for canonical pathways and GO terms in PDO versus tumour, PDX versus tumour and PDO versus PDX. Green and purple boxes indicate log₁₀(P-values) for down- or up-regulated gene sets, respectively. Tumour samples with a tumour purity ≥40% are shown. (c) Violin plots showing the expression of SCD and UGT1A10 in the OT cohorts (patient tumour purity ≥40%, PDX and PDO) and CCLE 2D colorectal cancer cell lines (see methods). Lines mark the 25%, 50% and 75% quantiles.

Figure 6. Overview of drug response in PDX and PDO models.

(a) Box plot showing the maximum inhibition (E_max) values for 14 drugs and for the staurosporine control in 35 PDO models. Irinotecan was tested on 21 PDOs. Range of drug efficiency (E_max>50%) or inactivity (E_max<50%) is indicated in blue and grey, respectively. (b) Box plot showing the potency of 14 drugs and of the staurosporine control based on log₁₀(IC₅₀) in 35 PDO models. The four response categories are coloured as indicated in the caption. (c) Hierarchical clustering representing the correlation of drug sensitivity patterns in PDOs based on the log₁₀(IC₅₀) values. (d) Box plot showing the sensitivity of 57 PDX models to 11 drugs based on tumour/control volumes (T/C in %) values: response categories are coloured as indicated in the caption. (e) Hierarchical clustering representing the correlation of drug sensitivity patterns based on T/C values in PDXs. (f) Waterfall plot showing cetuximab sensitivity for 53 PDX and indicating mutations and fusions in RAS/RAF and RTK pathways. Both T/C and RECIST results are indicated; models derived from the same patient are labelled with the same symbol above each sample name. The vertical dash line shows the cut-off between responders and non-responders. (g) Analysis of the activating ALK gene fusion in 353_T_XEN; upper part: schematics displaying chromosomal location, fusion partners, exon structure (grey and red boxes), fusion breakpoint (red line), validation by Sanger sequencing and RNA read coverage. Lower part: treatment response of OT PDX 353_T_XEN to crizotinib, cetuximab, AZD8931 and afatinib: mean tumour volume of treated and corresponding controls measured between day 8 and day 36 after treatment.

Figure 7. Bubble plot comparison of drug sensitivity profiles in 19 sibling pairs of PDX and PDO models.

Drug response categories are displayed on the axes (x=PDX, y=PDO). The size of each bubble is proportional to the number of model pairs in a given response category, where this number is indicated in the bubbles. The colours represent the eight drugs tested in both model systems, as indicated in the caption. The background colour in each square indicates the degree of concordance between the two models where grey indicate discordant drug response.

Figure 8. Molecular classification of response to 5-FU in PDX models.

(a) Hierarchical clustering of PDX samples based on the 5-FU response gene signature (n=251 genes). The drug sensitivity is given by T/C continuous values coloured as indicated in the caption (dark green: strong response, red: resistance). MSI samples, OT-NMF-groups and CMS classifications are coloured as indicated in the caption. (b) Hierarchical clustering of PDX samples based on the mini-classifier (n=14 genes) obtained by SVM-based machine learning on 49 PDX. (see Supplementary Data 14 for cross-validation). (c) Distribution of the predicted responders and resistant patient tumours (tumour purity ≥40%) in the three molecular groups (ECM/EMT, ASCL2/MYC, Entero/Goblets).

Figure 9. Molecular classification of response to EGFR inhibitors in PDX and PDO models.

Drug sensitivity is represented by T/C (PDX) or IC₅₀ (PDO) continuous values coloured as indicated in the caption (dark green: strong response, red: resistance). MSI samples, OT-NMF-groups and CMS classifications (SSP: upper row, RF: bottom row) and mutation status are coloured as indicated in the caption. (a,b) Hierarchical clustering of 33 PDO samples based on the drug response gene signatures for afatinib (a) or AZD8931 (b). (c,d,e) Hierarchical clustering of 49 PDX samples based on the drug response gene signatures for afatinib (c), AZD8931 (d) or cetuximab (e). Copy numbers (CN) of MYC (8q24.21) and AURKA (20q13.2) loci in PDX samples are shown below the cetuximab signature. The bars are colour-coded according to the CNV groups in the patients (Hypoploid, Hyperploid and MSI). Violin plots show the AURKA and MYC expression in responders versus resistant PDXs (response T/C≤25%, resistance T/C>25%). The expression values are shown as log₂ RPKM. Lines indicate the 25%, 50% and 75% quantiles, respectively.

Figure 10. Validation of the mini-classifier of cetuximab response in PDX.

The drug sensitivity is given by T/C continuous values coloured as indicated in the caption (dark green: strong response, red: resistance). MSI samples, OT-NMF-groups, CMS classification and mutations are coloured as indicated in the caption. (a) Hierarchical clustering of 48 PDX samples based on the mini-classifier (n=16 genes) obtained by SVM-based machine learning. (b) Principal component analysis (PCA) showing the classifier's ability to separate responders from non-responders. (c) Respective performances of the OT mini-classifier and of the RAS/RAF mutation status in predicting cetuximab sensitivity in the following cohorts OT-PDX (OT) (cross-validation), EPO PDX, the Gao et al., PDX (NV) and Khambata-Ford et al., primary tumours (KF). The mutation status was defined by mutations in codon 12 and 13 of KRAS or detected activating mutations in KRAS, BRAF or NRAS (BRAF mutations: V600E; KRAS/NRAS mutations: G12, G13, Q22, Q61, A146). In an additional setup, stable disease (SD) samples were excluded. The dot plot shows the specificity, sensitivity and balanced accuracy of the predictions. The 0.8 value is highlighted.with a black dotted line. (d) Same as in (c) but taking as validation cohort 164 samples merged from OT, EPO, NV and KF. (e) Distribution of the predicted responders and resistant patient tumours (tumour purity ≥40%) in the three molecular groups (ECM/EMT, ASCL2/MYC, Entero/Goblets).