An integrated omics approach highlights how epigenetic events can explain and predict response to neoadjuvant chemotherapy and bevacizumab in breast cancer

Abstract

Treatment with the anti-angiogenic drug bevacizumab in addition to chemotherapy has shown efficacy for breast cancer in some clinical trials, but better biomarkers are needed to optimally select patients for treatment. Here, we present an omics approach where DNA methylation profiles are integrated with gene expression and results from proteomic data in breast cancer patients to predict response to therapy and pinpoint response-related epigenetic events. Fresh-frozen tumor biopsies taken before, during, and after treatment from human epidermal growth factor receptor 2 negative non-metastatic patients receiving neoadjuvant chemotherapy with or without bevacizumab were subjected to molecular profiling. Here, we report that DNA methylation at enhancer CpGs related to cell cycle regulation can predict response to chemotherapy and bevacizumab for the estrogen receptor positive subset of patients (AUC = 0.874), and we validated this observation in an independent patient cohort with a similar treatment regimen (AUC = 0.762). Combining the DNA methylation scores with the scores from a previously published protein signature resulted in a slight increase in the prediction performance (AUC = 0.784). We also show that tumors receiving the combination treatment underwent more extensive epigenetic alterations. Finally, we performed an integrative expression-methylation quantitative trait loci analysis on alterations in DNA methylation and gene expression levels, showing that the epigenetic alterations that occur during treatment are different between responders and non-responders and that these differences may be explained by the proliferation-epithelial-to-mesenchymal transition axis through the activity of grainyhead like transcription factor 2. Using tumor purity computed from copy number data, we developed a method for estimating cancer cell-specific methylation to confirm that the association to response reflects DNA methylation in cancer cells. Taken together, these results support the potential for clinical benefit of the addition of bevacizumab to chemotherapy when administered to the correct patients.

Keywords: DNA methylation; bevacizumab; breast cancer; chemotherapy; epigenetics; multiomics.

Fig. 1

Flow chart of analyses indicating key results. The analyses were split in two main parts: prediction of response and epigenetic alterations during treatment. AUC, area under the curve; ChIA‐PET, chromatin Interaction Analysis by Paired‐End Tag sequencing; GSEA, gene set enrichment analysis; pCR, pathological complete response; TF, transcription factor.

Fig. 2

Identification of predictive signatures for response to neoadjuvant treatment with chemotherapy plus bevacizumab in ER‐positive breast cancer. (A) Receiver operating characteristic (ROC) curve showing specificity vs sensitivity for the leave‐one‐out cross‐validation probabilities for response. (B) Boxplots of DNA methylation levels in pre‐treated samples for example CpGs in the identified predictive signature. P‐value is determined using a two‐sided t‐test. (C, D) ROC curves showing specificity vs sensitivity in the PROMIX validation cohort (C) the “translated” gene expression signature, and (D) the hybrid protein‐methylation score. AUC, area under the curve; Se, sensitivity; Sp, specificity.

Fig. 3

Identification of biclusters of CpGs and genes altered during treatment with chemotherapy +/− bevacizumab. (A) Spectral co‐clustering of correlation coefficients of delta expression‐methylation quantitative trait loci (emQTL) using k‐means identifies two biclusters; heatmap shows ordered rows and columns, where blue represents negative correlation and red represents positive correlation. (B) Gene set enrichment analysis (GSEA) of genes in the two biclusters. (C, D) functional epigenomic enrichment (ChromHMM) analysis of CpGs (C), and transcription factor binding region (UniBind) enrichment analysis (D) of CpGs in Bicluster 1 and Bicluster 2 (upper and lower panels, respectively). P‐values (x‐axis) are calculated using hypergeometric tests and fill color (C, D) denotes fold enrichment above what is expected by chance. (E) Fold enrichment of overlap between CpG‐gene pairs in biclusters and chromatin loops determined by chromatin interaction analysis by paired‐end tag sequencing (ChIA‐PET) in MCF7; statistical significance (Bicluster 1: P = 8.6e‐30(*); Bicluster 2: P = 1) was determined by hypergeometric test of observed loops compared to possible loops within each biclusters.

Fig. 4

Delta DNA methylation of the identified biclusters after treatment with chemotherapy +/− bevacizumab and association to treatment response. Hierarchical clustering and heatmap of delta methylation values in Bicluster 1 (A) and Bicluster 2 (B); red is gain of methylation and blue is loss of methylation. Three patient clusters were in identified in both dendrograms. Patients (columns) are annotated with PAM50 subtype, estrogen receptor (ER) status, administration of bevacizumab and whether the patient achieved pathological complete response (pCR). (C–D, F–G) Boxplot showing average delta methylation of Bicluster 1 and 2 CpGs, respectively, plotted against achievement of pCR or residual cancer burden (RCB). Black dashed line denotes no change in methylation. Statistical significance is calculated using a t‐test. (E, H) Correlation between average delta methylation of Bicluster 1 and 2 CpGs, respectively, and the fraction of tumor left for all patients (black), combination arm (orange) and chemo arm (brown). Dashed lines denote linear regressions between average delta methylation and response. Statistical significance is calculated using Pearson correlation test.

Fig. 5

Delta gene expression of the identified biclusters after treatment with chemotherapy +/− bevacizumab and association to treatment response. (A, B) Hierarchical clustering and heatmap of delta expression values; red is gain of expression and blue is loss of expression. Three patient clusters were in identified in both dendrograms. Patients (columns) are annotated with PAM50 subtype, estrogen receptor (ER) status, administration of bevacizumab and whether the patient achieved pathological complete response (pCR). (C–D, F–G) Boxplot showing average delta expression of Bicluster 1 and 2 genes, respectively, plotted against achievement of pCR or residual cancer burden (RCB). Black dashed line denotes no change in expression. Statistical significance is calculated using a t‐test. (E, H) Correlation between average delta expression of Bicluster 1 and 2 genes, respectively, and fraction of tumor left for all patients (black), combination arm (orange) and chemo arm (brown). Dashed lines denote linear regressions between average delta expression and response. Statistical significance is calculated using Pearson correlation test.