BRAF V600E Mutant Colorectal Cancer Subtypes Based on Gene Expression

Colorectal cancer is the third leading cause of cancer-related death in the United States (1). The BRAF gene is activated by mutation in 50% of melanomas and 10% of colorectal cancer (2). In the MAPK pathway, KRAS activates BRAF, which then activates MEK kinase (MEK1 and MEK2) and finally ERK (3). However, the nonrandom occurrence of BRAF mutations in certain disease subtypes such as hypermethylated right-sided colorectal cancer or specific thyroid carcinoma subtypes suggests that additional tumor features and alterations are associated to the presence of V600 BRAF and will determine the final signal output (2). Moreover, the gene expression patterns of KRAS mutant and BRAF mutant colorectal cancer are typically very different from each other (4, 5). Of interest, the oncogenic contribution of mutated BRAF may vary between tumor types, as suggested by the very heterogeneous clinical benefit seen for BRAF inhibition strategies between melanoma and colon (6, 7). This heterogeneity in drug resistance indicates cancer-type specificity of the biology underlying the BRAF mutation. In addition, the recent outcome data of trials targeting selected BRAF V600 mutated colorectal cancer cohorts with monotherapy (8, 9), double-targeted therapy (8, 10–13), and triple-targeted therapeutics (13–15) have shown a surprisingly strong heterogeneity in response, suggesting that further biological subdivisions may exist even within BRAF V600 mutant colorectal cancer.

The BRAF V600E mutation (from now one referred to as BM) results in a constitutively active protein (2) and accounts for 96% of BRAF V600 mutations (16). More than 30 other oncogenic BRAF mutations were found. Although most of them also cause constitutive kinase activity, other mutations produce kinase-dead alterations that nevertheless can lead to oncogenic signals via CRAF activation (17, 18). The prognostic role of non-V600E BRAF mutations is still uncertain as they are not routinely tested for in clinical trials. The BRAF mutation is associated with other clinical factors such as the microsatellite instability (MSI) status (58% of BM patients are MSI), tumor site (81% of BM patients have right-sided tumors), and gender (62% of BM patients are females; refs. 19, 20).

In general, colorectal cancer patients have so far been treated irrespectively of clinical or molecular features (except exclusion of anti-EGFR therapy for RAS-activated tumors) and display heterogeneous treatment response. Major efforts have been devoted to defining molecular colorectal cancer subtypes with the increasing use of transcriptomics data (21–28), to which we have extensively contributed (22, 29). Colorectal cancer classification is a highly active field of research as it aims to stratify patients for therapeutic interventions and provides more precise assessment of heterogeneous trial outcomes. Recently, a concerted effort by the Colorectal Cancer Subtyping Consortium (CRCSC) succeeded in unifying colorectal cancer subtyping (29). This led to subdivision of colorectal cancer tumors into four consensus molecular subtypes (CMS), capturing the main overall gene expression variability in colorectal cancer (29). The CRCSC system aimed at consolidating the variation that had been consistently modeled in the previously subtype publications, without investigating evidence for additional variation. The vast majority (> 70%) of BRAF mutants were classified into the same subtype (CMS1) whereas 17% fell into CMS4 and 6% into CMS3. We decided to perform a dedicated analysis to better characterize the heterogeneity within BRAF mutants and collected a large set of 218 BM patients with gene expression, basic mutation, and outcome data for further study. Here we report that BM patients can be separated into two different gene expression subtypes with distinct molecular patterns and potential different therapeutic targets.

The database used in this study was constituted by combining the following independent cohorts: colorectal adenocarcinomas from The Cancer Genome Atlas (TCGA; ref. 21), the Pan-European Trials in Alimentary Tract Cancers (PETACC-3; refs. 5, 30), three different datasets from Agendia (GSE42284; ref. 31), ICO (4, 29), and VHB (4, 29), Marisa and colleagues (GSE39582; ref. 24), Schlicker and colleagues (GSE35896; ref. 28), and two datasets from the Ludwig Institute for Cancer Research in Melbourne, GSE75315 and GSE75316. The number and origin of patients with some of their clinicopathologic characteristics are highlighted in Figure 1. Detection of clinicopathologic features, MSI status, KRAS, and BRAF mutation status was performed as described in the original reports describing the cohorts (4, 5, 21, 24, 28, 30, 31). The mutational analysis collected by each study is summarized in Supplementary Table S1. Briefly, KRAS mutants include mutations at residue G12, G13, and Q61; BRAF mutants only include the V600E mutation and PIK3CA mutants includes all known activating mutations.

Unless otherwise mentioned, all analyses in this study were performed using R software (version 3.2.2). A reproducible-research script and data can be found at our institutional website (32). Level 3 TCGA RNA-sequencing expression data were obtained from the TCGA portal. GSE39582 and GSE35896 were obtained preprocessed from the Gene Expression Omnibus (GEO) website. PETACC-3 expression data had been processed as described earlier (5). Probe levels expression datasets were processed as follows: for probes mapping to a unique gene, only the probe displaying the largest variation was conserved. Gene expression data from all cohorts were merged at the gene level, which resulted in a dataset of 9,328 genes for 218 BRAF V600E mutant patients. The ComBat function from the sva package was used to normalize and correct for batch effects across datasets (33). The similarity structure of the tumor profiles was then assessed by principal component analysis (PCA) and hierarchical clustering.

Cluster analyses were performed using correlation distance metrics and the average linkage agglomeration algorithm (R package nclust version 1.9.0; ref. 34). Non-negative matrix factorization (NMF) was done with the NMF package (version 0.20.5) and standard strategies. Detailed descriptions of the NMF-based subtyping is described in the Supplementary Methods. Statistical analyses assessing the association between BM subtypes and clinicopathologic factors were performed by Fisher exact tests. P-values were corrected for multiple comparisons using the Bonferroni method. The MSI-adjusted clustering analysis was performed as described in the Supplementary Methods.

We constructed a classifier for the two BM groups by training a penalized logistic model with the subtype labels given by the selected two subgroups NMF solution, using the glmnet package (version 2.0.2). We obtained a classifier with nonzero coefficients for 44 genes listed in the Supplementary Table S2. Methodological details are described in the Supplementary methods. To determine genes that characterizes the difference between BM1 and BM2 without being confounded by relevant clinical factors (MSI status, site, and gender), which were not perfectly balanced between BM subtypes, we fitted a multivariable linear regression including the three potential confounders as additional explanatory variables in the model (35). To extract the genes specifically associated to the BM subtypes, we selected the genes with a P-value inferior to 5 × 10⁻⁸ [corresponding to a Benjamini–Hochberg FDR (BH-FDR) of 0.002] for the effect of the BM subtype variable on gene expression, which led to the identification of 476 genes.

Pathway analyses were performed by Gene Set Enrichment Analysis (GSEA) and are described in the Supplementary Methods.

Analyses performed on external public data for DNA methylation, proteomics, and cell line gene expression are described in the Supplementary Methods, along with methods for survival analyses. The public cancer cell lines data were from refs. 36 and 37.

We built a gene expression dataset of 218 BRAF V600E mutant (BM) colorectal cancer patients from different clinical trials (Fig. 1) to study gene expression profiles in this relatively rare subgroup. Although the profiles were obtained using different technologies and using fresh or FFPE specimen as source material, they did not show any evident cohort-bias clustering by PCA (Supplementary Fig. S1A) or hierarchical clustering (Supplementary Fig. S1B).

To assess if the tumor data presented a cluster structure, we applied a class-discovery approach based on NMF. The NMF algorithm is stochastic; we computed 200 solutions for models with a factorization rank of two to six components. We judged the quality with two classic criteria: the cophenetic correlation coefficient (CCC, a measure of the overall stability of sample clustering over multiple solutions with a given rank) and the silhouette width (SW; a measure of the similarity of each sample to the samples in the same cluster relative to the similarity to the other samples; ref. 38). Very informative is also visual inspection of the consensus matrix (Fig. 2A), which suggests that the tumor set falls naturally into two separate groups. In agreement with this interpretation, the CCC and the average SW decrease significantly with increasing rank, and particularly from rank two to rank three. For rank two the CCC and SW were close to their maximum value of one, indicating, on an absolute scale, an appropriate clustering into two subtypes of tumors, whereas these criteria are not defined for rank one (the original data). A gene-randomized control dataset was not split at all (Fig. 2A). We will refer to these subtypes as BM1 (BRAF mutant 1) and BM2 (BRAF mutant 2). Of 218 BM colorectal cancer, 69 were classified into BM1 and 149 into BM2 (BM1:BM2 ratio of approximately 1:2). An algorithm that performs an unsupervised biclustering of genes and tumors confirmed that the two BM subtypes cluster separately (Supplementary Fig. S2A). These results indicated that the colorectal cancer BRAF V600E mutant population is constituted by two main subtypes, characterized by specific gene expression patterns.

We then studied if these subtypes correspond to groups defined by known clinical factors. Fisher exact test revealed a borderline statistically significant association between BM subtype and MSI status (P = 0.04) and a trend toward association with gender (P = 0.11; Fig. 2B-C; Supplementary Table S3). BM2 was slightly enriched in MSI (65% of MSI in BM2 vs. 43% in BM1) and female patients (68% of female in BM2 vs. 51% in BM1; Fig. 2B and C), but these differences do not drive the subtype split, as shown below in the conditional model. Activating PIK3CA mutations, primary tumor site and stage were homogeneously distributed across BM subtypes (Fig. 2B). The BM1:BM2 ratio of approximately 1:2 was observed across all cohorts, indicating that the NMF clustering was not reflecting batch or cohort effects (Fig. 2D).

The observed association between BM subtypes and MSI status raised the question whether the clustering itself or many of the gene expression differences between BM1 and BM2 were due to the different proportion of MSI tumors. To answer this question, we performed two additional analyses. First, we did an NMF clustering after adjusting the expression data for effects associated to MSI status (see Materials and Methods). The subtypes obtained were almost identical [Cohen κ value of 0.82 (95% CI = 0.73–0.90)], indicating that the particular split into two groups is not driven by differences between MSI and MSS tumors (Supplementary Fig. S2B). Second, we fitted a multivariable linear regression for all genes using BM subtype, gender, MSI status, and tumor location as explanatory variables (Fig. 2E; Supplementary Table S4). We found that 476 genes were differentially expressed between BM1 and BM2 and were specific to the BM subtype variable (Fig. 2E), 48 genes were strongly associated with MSI status and 8 with gender. Importantly, only one gene was associated with both MSI status and BM subtype, indicating the independence of these two variables at the gene expression level.

The CRCSC subtype consensus consists of four subtypes called CMS1 (characterized by MSI and immune patterns), CMS2 (chromosomal instability and WNT activation), CMS3 (metabolic pattern), and CMS4 (epithelial–mesenchymal transition, EMT). In this study, the vast majority of BRAF mutant patients (70%) were classified into CMS1 whereas only a few were found in CMS2 (2%), CMS3 (5%), and CMS4 (17%). More than 50% of CMS1 tumors are non-BRAF mutants, indicating that by focusing our analysis on the BRAF mutant population, we might capture more variation specifically associated to the BRAF mutation. Interestingly, all CMS4 BRAF mutants are classified as BM1, whereas CMS1 BRAF mutants are distributed into both BM1 and BM2 (Fig. 2B and C). These results are not in contradiction with the CMS classification but instead indicate that our classification is refined for BM patients by capturing additional transcriptomics variation within the BRAF mutant population that is part of CMS1. Altogether, these observations indicate that the BM subtypes are not related to an already known classification system and thus are new yet uncharacterized groups.

We conducted gene set enrichment analyses (GSEA) to characterize the molecular pathways specific for the BRAF mutant subtypes. Using the 476 subtype-specific genes, we scored signatures of hallmark processes collected in the MSigDB molecular signature database (39) using the single sample GSEA (ssGSEA) method that indicates the degree to which the genes in a particular gene set are coordinately up- or downregulated within a sample (40). Prominent signatures scoring high in BM1 were EMT-related processes (EMT, apical junction, and myogenesis), KRAS signaling, and immune response (Fig. 3A; Supplementary Fig. S3). High scores in BM2 were achieved by cell-cycle and cycle checkpoints-related processes such as target genes of E2F transcription factors and genes involved in the G₂ to M (G₂–M) transition (Fig. 3A; Supplementary Fig. S3). Figure 3A indicates that several molecular signatures were characteristic for BM1, respectively for BM2, in a direct comparison between these two groups but not necessarily in comparison to other colorectal cancer tumors. To broaden the analysis to include all colorectal cancer, we compares ssGSEA signature scores between four groups: BRAF mutant subtype, KRAS/BRAF double wild-type (WT2), and KRAS mutant colorectal cancer (see Materials and Methods). The KRAS signature was not only enriched in BM1 compared to BM2 but also compared with KRAS mutant colorectal cancer (Fig. 3B; Supplementary Fig. S4). BM2 displays similar levels of KRAS signaling than WT2 colorectal cancer. Cell cycle–related processes and glycolysis were not only overactivated in BM2 compared with WT2 and KRAS mutant colorectal cancer but also repressed in BM1 (Fig. 3B; Supplementary Fig. S4). Similarly, EMT-related processes were enriched in BM1 but also repressed in BM2 (Fig. 3B; Supplementary Fig. S4).

With the aim of characterizing the BM subtypes in their globality, we also conducted pathway analyses using the whole gene expression profile rather than by focusing on the 476 differentially expressed genes. In addition to the pathways found with the reduced list of genes, BM1 displays an overall stronger immune profile emphasized by activation of pathways such as IL2/STAT5, TNFα signaling via NF-κB, IL6/JAK/STAT3, and allograft rejection (Supplementary Fig. S5). Angiogenesis and TGFβ processes were also enriched in BM1 compared to BM2. BM2 enrichment in metabolic process was more evident by looking at overall gene expression. This full transcriptomics analysis also revealed that both BM groups display less WNT activation than WT2 and KRAS mutant colorectal cancer, and that BM2 has less WNT activity than BM1 (using the M5895 signature). Of note, this difference in WNT activation was not present by using the 476 BM specific genes (Supplementary Fig. S5).

To confirm the consistency of the signature enrichment, we performed ssGSEA with manually selected signatures covering the biological processes found to be enriched in Fig. 3A (Fig. 3C). The vast majority of KRAS and EMT-related signatures suggested activation of these processes in BM1. The same was true for cell-cycle–related signatures that were scoring high in BM2. Noteworthy, there was minimal gene overlap across this set of gene signatures (Supplementary Fig. S6).

The mTORC1 signature (MSigDB ref: M5924) was scoring high in BM2. An mTORC1 complex activation can be the consequence of the activity of several other pathways including EGFR, MAPK, AKT, cytokines, hypoxia, and WNT (41). In turn, mTORC1 complex activation can lead to 4EBP1 and S6K pathways activation (41). Thus, to understand better the BM2 enrichment in mTORC1 signature, we have analyzed in greater details the scoring of some of these pathways in the BM subtypes. Surprisingly, we found that mTOR, 4EBP1, and two WNT signatures were enriched in BM1 instead of BM2 (Fig. 3C). To elucidate this discrepancy, we checked the overlap of these signatures and found none (Supplementary Fig. S6). We therefore checked for KEGG biological process enrichment in these two signatures and found that M7561 (52 genes) was indeed highly representative of the mTOR pathway whereas M5924 (200 genes) was not and was instead rather enriched for genes related to cell-cycle and metabolic processes (Supplementary Table S5). The pathway analyses, therefore, indicate that BM1 is characterized by mTOR and RAS signaling, whereas BM2 by cell-cycle and metabolic processes. The WNT signaling enrichment found in BM1 was not consistent across all signatures and should be interpreted cautiously.

To study in greater details the immune profile of BM subtypes, we used the recently reported tool called CIBERSORT that estimates the immune landscape based on gene expression profiles. This analysis revealed that BM1 tumors were significantly more strongly infiltrated by macrophages and monocytes and significantly less by resting dendritic cells whereas BM2 displayed more mast cell content (Supplementary Fig. S7; ref. 42).

Transcriptomics data do not account for posttranslational modifications. Proteomics reverse phase protein array (RPPA) and phosphoproteomics data are publicly available for some colorectal cancer of the TCGA, including 6 BM1 and 21 BM2 TCGA patients. The results of a differential expression analysis were consistent with the reported pathway analyses (Fig. 4A). In particular, AKT Serine 473 and 4EBP1 Threonine 70 were significantly more phosphorylated in BM1, whereas eukaryotic translation initiation factor 4E (eIF4E), another player of the mTOR pathway, was downregulated in BM1 (Fig. 4A and B; linear regression P-value < 0.05). Differentially expressed proteins also included cell cycle–related proteins such as CDK1, cyclin D1, and ATM, which confirms our observation that BM2 displays cell-cycle misregulation (Fig. 4A and B). Noteworthy, other proteins were associated with BM subtypes such as the BH3-only proapoptotic BIM protein, in line with the high-scoring apoptosis signatures seen in BM1 (Supplementary Fig. S4–S5). BM1 also displayed an enhanced inflammatory response, which is corroborated by the differential expression of proteins such as SYK that transmits signals in B and T cells and STAT5α whose expression correlates with immune activation.

The association between BRAF mutation and the CpG island methylation phenotype (CIMP) is well documented (43). In TCGA methylation data, 4% of the probes were significantly more methylated in BRAF mutant than in wild-type colorectal cancer (Supplementary Fig. S8A); however, we found no differences between the global methylation patterns of BM1 and BM2 (Supplementary Fig. S8B). In view of the overrepresentation of MSI colorectal cancer in BM2, which might be associated with the MLH1 hypermethylation (44), we studied the methylation status of MLH1 in TCGA data (Supplementary Fig. S8C). MLH1 hypermethylation was visually obviously associated to MSI status as expected, but not to BM subtypes. The few differentially methylated probes between BM1 and BM2 (Supplementary Fig. S8D) were not obviously linked to the pathway analyses reported above.

We compared overall and relapse-free survival between BM subtypes using the Cox proportional hazard model. We observed a trend for the patient survival in the BM1 group to be poorer than that in the BM2 group (Fig. 5 and univariable test in Table 1; see patient characteristics in Supplementary Fig. S9). Multivariable testing suggested that MSI status (MSI and MSS) is the dominant prognostic factor whereas BM information does not add prognostic power in a model that includes MSI status (multivariable test in Table 1) or the remaining difference is too small for the power of this analysis to be statistically significant (BM subtype + MSI status: N = 140; events = 42; BM subtype + stage: N = 166; events = 50). Subset survival analysis in MSI status and in tumor-stage subgroups did not show a prognostic effect for BM (subset analysis test in Table 1).

Therefore, the data offer no statistically significant evidence that the survival in BM1 and BM2 is different and the observed trend might be due to the unbalanced representation of MSI and MSS patients (Fig. 5 and Table 1).

The Cancer Cell Line Encyclopedia (CCLE) and the Sanger institute and GlaxoSmithKline (GSK) have published gene expression profiles and drug screening results for a wide panel of cancer cell lines. We used these resources to investigate whether BM subtypes might differ in drug-response patterns. We built a penalized logistic model classifier for BM1/BM2, with a penalty coefficient that minimizes the misclassification error estimated by cross-validation. The classifier uses 44 genes (Supplementary Table S2) and has an estimated misclassification error of 11% (Supplementary Fig. S10). We applied this classifier to colorectal cancer cell lines: SW1417 and C2BBe1 were classified as BM1 whereas HCT15, LS1034, SW1116, RKO, SW48, and SW837 were classified as BM2, consistently in both datasets (supplementary Table S6). Twenty-five percent of these cell lines harbor the BRAF mutation and no evident enrichment toward BM subtypes was observed (one out of two BM1 cell lines and one out of six BM2 cell lines harbor the BRAF mutation; Supplementary Fig. S11).

Next, we analyzed how these cell lines respond to selected drugs. According to Sanger data, a differential analysis revealed that BM1 cell lines are more sensitive to inhibition of BRAF (PLX4720; linear regression unadjusted P-value: 0.075), BCL2 (ABT-263; linear regression unadjusted P-value: 0.011), and MEK (CI-1040; linear regression unadjusted P-value: 0.101) compared to BM2 (Supplementary Table S7 and Supplementary Fig. S11). However, BM2 cell lines are more sensitive to CDK1 inhibition (by RO-3306) compared to BM1, although not statistically significantly. The analysis of CCLE data showed that an inhibitor of SRC (AZD0530), a protein that is upstream of RAS, RAF, PI3K, AKT, and MEK, is more efficient in killing BM1 cells compared to BM2 (linear regression unadjusted P-value: 0.076). The HSP90 chaperone 17-AAG displayed modest activity toward BM2 cells in both Sanger and CCLE data. These results support the notion that different pathways are preferentially activated in the two BM subtypes and allows prospective testing of novel drug combinations in stratified sets of models and patients.

The CRCSC recently published a consensus classification system that unifies previously reported colorectal cancer subtypes based upon gene expression profiles (29). It consists of four subtypes called CMS1 to CMS4. In the largest yet reported colorectal cancer BRAF mutant patient cohort we describe here, we found that the vast majority classifies as CMS1 and only a few in any other groups. This does not explain the heterogeneity seen among BRAF mutant patients at the clinical level. The gene expression profiles segregate these carcinomas naturally into two subtypes, BM1 and BM2, which differ by the expression of a group of 476 genes, independently from factors such as tumor site, MSI, and gender. The two subtypes displayed difference in overall and relapse-free survival, although not significantly. MSI is well known to confer good prognosis independently of other clinical factors (45). One explanation for the survival difference between subtypes resides in the enrichment of BM2 with MSI patients as shown by a Cox multivariable analysis. However, although MSI status is a confounding factor for survival, our data indicate that it does not drive the difference in gene expression. Enrichment of BM1 in EMT signatures and MSS patients, which are associated with poor prognosis, correlates with poor survival of patients with a CMS4 consensus subtype colorectal cancer (29). Our survival analysis has some limitations: there are slight intercohort differences in the definition of events in relapse-free survival, which include secondary primary tumors in one cohort, in the computation of time to event, and in the distribution of censoring times (Supplementary Fig. S9). Moreover, due to the low number of patients and events in subset analyses, the statistical power of these analyses is relatively low.

Detailed pathway analysis of the 476 differentially expressed genes and of the overall transcriptome as well revealed that the BM1 subtype is highly active in KRAS/mTOR/AKT/4EBP1 signaling and in genes associated with macrophage infiltration and EMT and the BM2 subtype in cell cycle and cycle checkpoint associated genes. Proteomics studies have revealed high levels of phosphorylation for 4EBP1 threonine 70 and AKT serine 473, which is in line with a high RAS signature. As a key effector of AKT-dependent oncogenic activation in tumors, 4EBP1 is released from its binding to eIF4E when phosphorylated by AKT, which subsequently initiates protein translation (46). Not surprisingly, phosphorylation of 4EBP1 at the threonine 70 residue is associated with poor survival in melanoma (47).

BRAF inhibition is known to initiate a feedback loop via reactivation of EGFR (48). As a result, single therapy with BRAF inhibitors has not been effective in colorectal cancer. An approach combining BRAF inhibition with EGFR blockade has been tested and showed promising results (48). Interestingly, BM1 subtype is enriched for an EGFR signature, derived from stable overexpression of ligand activated EGFR (M2634 in Fig. 3C). Nevertheless, we have found EGFR expression levels did not differ significantly between BM1 and BM2 (data not shown). Thus, this makes it unlikely that the pattern of resistance to BRAF inhibition differs between the two BM subtypes.

The CCLE and the Sanger institute have published gene expression profiles of a large panel of cancer cell lines. In addition, they screened these cell lines for drug sensitivity. We used these resources in an attempt to predict drug response according to BM subtype. Unfortunately, the very low number of cell lines classified as BM1 did not allow sufficient statistical power and therefore this analysis remains only exploratory. Nonetheless, our data suggest that BM1 and BM2 patients might benefit from different targeted approaches. A clinical trial assessing the combination of BRAF inhibitors with EGFR inhibitors and PI3K inhibitors (NCT01719380 in clinicaltrials.gov) is currently ongoing and its results are very much awaited. Our data suggest that BM1 patients may benefit from this combination (Fig. 4; Supplementary Fig. S11). An obvious link between BRAF activation and cell-cycle activation is the MAPK pathway (49). Another clinical trial evaluating the combination of BRAF inhibitors with EGFR inhibitors and MEK inhibitors (NCT01750918 in clinicaltrials.gov) is currently in progress with equally eagerly awaited results. Our data suggest that this combination might be effective in patients with a BM2 subtype colorectal cancer (Fig. 4; Supplementary Fig. S11).

Because BRAF mutant colorectal cancer cell lines and tumors do not respond efficiently to current anti-BRAF treatments such as Vemurafenib and PLX4720, strategies to combine with other drugs have recently been proposed (48, 50). Among them, combination of PLX4720 with the AZD6244, a MEK1/2 inhibitor and MK-2206, an AKT inhibitor, or with AZD7762, a CHK1/2 checkpoints inhibitors appeared highly effective in colorectal cancer cell lines (50). In this respect, the finding that BM1 is characterized by signatures related to RAS, a player of the MAP kinase signaling upstream of MEK, as well as related to AKT and 4EBP1, holds promise for treatment of patients with a BM1 tumor with combination of anti-BRAF, anti-MEK, and anti-AKT/mTOR/4EBP1 drugs. The combination of PLX4720 with the checkpoint kinase 1/2 (CHK1/2) inhibitor called AZD7762 strongly induced cell death in several BRAF mutant cell lines (50). This suggests that drugs targeting cell-cycle components, such as CHK1/2 and Aurora kinases, might be effective in BM2 patients. These hypotheses, however, need to be further explored.

Understanding what drives the difference between BM1 and BM2 is challenging. Because both subtypes harbor the BRAF mutation, it cannot itself be determining their differences. The analysis of the TCGA cohort did not reveal any differences in mutational events between BM1 and BM2 (data not shown). However, this analysis lacks statistical power due to low number of patients. Thus, it is difficult to say whether the nature of the “driving event” occurs at the mutational, epigenetic, or is a question of the state of the cell-of-origin of the cancer or of its interaction with the tumor stroma, for example.

We conclude that BM colorectal cancer can be segregated into two different subtypes based upon gene expression profiles. Further studies are needed to further elaborate the potential clinical relevance of these subtypes, notably in terms of different combinations of targeted therapies.

M. Delorenzi reports receiving a commercial research grant from Bayer and has ownership interest (including patents) in Novartis and Roche. No potential conflicts of interest were disclosed by the other authors.

Conception and design: D. Barras, I.M. Simon, A.D. Roth, S. Tejpar, M. Delorenzi

Development of methodology: D. Barras, J. Guinney, S. Tejpar, M. Delorenzi

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): O.M. Sieber, R.N. Jorissen, C. Love, P.L. Molloy, I.T. Jones, S. McLaughlin, P. Gibbs, I.M. Simon, S. Tejpar

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): D. Barras, E. Missiaglia, P. Wirapati, R.N. Jorissen, F.T. Bosman, S. Tejpar, M. Delorenzi

Writing, review, and/or revision of the manuscript: D. Barras, E. Missiaglia, O.M. Sieber, R.N. Jorissen, C. Love, P.L. Molloy, I.T. Jones, P. Gibbs, J. Guinney, I.M. Simon, A.D. Roth, F.T. Bosman, S. Tejpar, M. Delorenzi

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): D. Barras, O.M. Sieber, P. Gibbs, S. Tejpar

Study supervision: F.T. Bosman, S. Tejpar, M. Delorenzi

We thank Josep Tabernero, Ramon Salazar, and Gabriel Capella for collecting patient material (4, 26, 31). We thank all individuals involved in the PETACC-3 trial. We acknowledge the Victorian Cancer Biobank for the provision of specimens and BioGrid Australia for providing clinical data and Delphine Grun for testing the program code allowing figure reproducibility.

This work was supported by the Medic Foundation, the Ludwig Institute for Cancer Research, the National Health and Medical Research Council of Australia (APP1050177, APP1062226), and the Victorian Government's Operational Infrastructure Support Program. Sabine Tejpar is supported by the KU Leuven GOA/12/2106 grant, the EU FP7 Coltheres grant, the Research Foundation Flanders-FWO, and the Belgian National Cancer Plan.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.