Evaluation of biological pathways involved in chemotherapy response in breast cancer

Abstract

Introduction

Our goal was to examine the association between biological pathways and response to chemotherapy in estrogen receptor (ER)-positive and ER-negative breast cancers separately.

Methods

Gene set enrichment analysis including 852 predefined gene sets, was applied to gene expression data from 51 ER-negative and 82 ER-positive breast cancers that were all treated with a preoperative paclitaxel, 5-fluoruracil, doxorubicin and cyclophosphamide chemotherapy.

Results

Twenty-seven (53%) ER-negative and 7 (9%) ER-positive patients had pathologic complete response (pCR) to therapy. Among the ER-negative cancers, a proliferation gene signature (FDR q = 0.1), the Genomic Grade Index (FDR q = 0.044) and the E2F3 pathway (FDR q = 0.22, p = 0.07) signature were enriched in the pCR group. Among the ER-positive cancers, the proliferation signature (FDR q = 0.001) and the Genomic Grade Index (FDR q = 0.015) were also significantly enriched in cases with pCR. Ki67 expression, as single gene marker of proliferation, did not provide the same information as the entire proliferation signature. An ER-associated gene set (FDR q = 0.03) and a mutant p53-gene signature (FDR q = 0.0019) were enriched in ER-positive cancers with residual cancer.

Conclusion

Proliferation- and genomic grade-related gene signatures are associated with chemotherapy sensitivity in both ER-negative and -positive breast cancers. Genes involved in the E2F3 pathway are associated with chemotherapy sensitivity among ER-negative cancers. The mutant p53-signature and expression of ER-related genes were associated with lower sensitivity to chemotherapy in ER-positive breast cancers only.

Introduction

Drug resistance is caused by multiple mechanisms that operate simultaneously in cancers. A large number of biological functions including transmembrane trafficking, DNA-repair, stress response, proliferation and apoptosis may affect the sensitivity of a cell to chemotherapy. Other, yet to be identified mechanisms may also play a role. Preoperative chemotherapy provides an attractive clinical setting to study mechanisms of drug resistance in patients.

Chemotherapy before surgery is used in the treatment of newly diagnosed, stage II-III breast cancers because it frequently reduces tumor size and improves surgical outcome [1]. Among patients who receive preoperative chemotherapy, up to 25-30% (depending on the type of treatment) experience complete eradication of the invasive cancer in the breast and regional lymph nodes after completion of 3-6 months of chemotherapy [2]. This favorable response is called pathologic complete response (pCR) and it indicates an extremely chemotherapy sensitive tumor and also heralds excellent long-term cancer-free survival [3]. We previously conducted a pharmacogenomic study that included 133 patients with newly diagnosed breast cancer who received preoperative chemotherapy with paclitaxel followed by 5-fluorouracil, doxorubicin and cyclophosphamide. All patients underwent a one-time, pre-treatment fine needle biopsy of the cancer for gene expression analysis. The goal of the study was to discover gene-expression based predictors of pCR. Our previous analysis focused on discovering the best possible multi-gene predictor without considering the function of any of the genes [4]. The goal of the current analysis is to examine an association between known biological pathways and response to chemotherapy.

Lists of genes (i.e. gene sets) that represent various biological pathways were assembled from the literature. We used gene set enrichment analysis (GSEA) to examine the correlation between these a priori defined gene sets and chemotherapy response [5]. Clinical experience as well as molecular analysis of breast cancers indicate that estrogen receptor (ER)-positive and ER-negative cancers are two different types of neoplastic diseases of the breast therefore [6, 7]. It is plausible that different molecular mechanisms may determine response or resistance to chemotherapy in these 2 types of breast cancers. Therefore, we performed our analysis separately for ER-positive and ER-negative cancers.

Materials and methods

Patients

This study included 51 ER-negative and 82 ER-positive tumors from patients with newly diagnosed stage I-III breast cancer. Each patient had a fine needle aspiration of the cancer before starting chemotherapy. These needle aspiration samples contain approximately 80% neoplastic cells and few or no stromal cells or normal breast epithelium [8]. All patients were treated with 6 months of preoperative chemotherapy with paclitaxel followed by 5-fluorouracil, doxorubicin and cyclophosphamide. Patients underwent surgery after completion of chemotherapy and the resection specimens were examined by a pathologist to measure residual cancer. For the purpose of our analysis, tumor response was dichotomized as pCR, defined as no residual invasive cancer in the breast and lymph nodes, or residual disease (RD) that included patients with any degree of invasive cancer that survived preoperative chemotherapy. The reason for this dichotomization was because pCR is a strong surrogate for long term cancer-free survival and therefore a marker of long-term benefit from therapy [2, 3]. It remains unknown to what extent patients who achieve less than pCR benefit from chemotherapy in terms of improved survival. This categorization of pathologic response allowed us to compare biological pathways between cancers with extreme chemotherapy sensitivity (pCR) and the rest (RD). There were not enough cases with tumor progression during treatment in our study to form a third group including extreme chemotherapy resistant tumors.

ER-status was determined from routine pathological assessment by immunohistochemistry. The cut-off for ER-positivity was ≥10% positive tumor cells following standard clinical practice. This study was approved by the institutional review boards (IRB) of MDACC and all patients signed an informed consent for voluntary participation. Clinical characteristics of the patients are presented in Table 1.

Table 1 Patient characteristics

Gene expression analysis

Gene expression profiling was performed by using Affymetrix U133A Gene Chips following standard operating procedures as described previously [4]. We normalized the gene expression data using dChip V1.3 software [9] to a single reference array. The normalized gene expression values were transformed to a log₁₀ scale for further analysis. The complete microarray data is available at the M.D. Anderson Cancer Center bioinformatics web site [10]. To identify differentially expressed genes between cases with pCR and RD we performed unequal variance t-test on each probe set. Because of the multiple comparisons many low p-values are expected by chance alone. Under the null hypothesis that no genes provide useful information, the distribution of p-values should be uniform. If, on the other hand, some genes do provide useful information about predicting response, we would expect an overabundance of small p-values (above what chance might produce). We can capture this situation by modeling the distribution of the p-values as a beta-uniform mixture (BUM). This analysis was used to estimate false discovery rates (FDR) that accompany particular p-values derived from t-test [11]. All analysis was performed using the R package (version 2.3.1].

Gene set enrichment analysis

Gene set enrichment analysis (GSEA) was applied to assess the association between pCR, RD, and 852 distinct a priori defined gene sets. The goal of the GSEA is to determine whether members of a particular gene set (i.e. a list of 15-500 probe sets that correspond to genes that define a biological pathway) tend to occur toward the top or the bottom of a rank ordered gene list including all gene expression measurements [5]. We ranked all probe sets based on their correlation with pCR. Three groups of gene sets were tested, the first included 319 distinct gene sets corresponding to probes associated with 295 different cryptogenic bands on 24 chromosomes, the second included 522 different gene sets corresponding to genes involved in various metabolic and signaling pathways. A detailed description of these gene sets and how they were assembled was presented by Subramanian et al [5]. The third group contained 11 gene sets of various oncogenic or drug resistance related pathways. These included 5 distinct oncogenic pathways that were defined as genes overexpressed in normal human mammary epithelial cells transfected with Myc, Ras, E2F3, β-catenin and Src oncogenes, respectively and were described by Bild et. al [12]. One mutated p53-associated gene set that was defined as genes overexpressed in human breast cancer with known p53 mutation and was previously described by Miller et al [13]. It has been suggested that this gene expression signature can distinguish cancers with wild type and mutant p53 and it may outperform direct p53 gene sequencing as a predictor of prognosis and therapeutic response. We also assessed one ER-associated gene set that contained genes that were most highly co-expressed with the ER gene in human breast cancer microarray data developed by Symmans et al [14]. This gene set did not include the ER gene itself. We also examined the genomic grade index (GGI) that represents genes that are differentially expressed between low grade and high grade human breast cancers and were identified by Sotiriou et al [15]. One prognostic signature that was derived by comparing gene expression profiles of tumors that recurred with those that did not was also tested. This prognostic signature was first reported by Wang et al [16]. A proliferation signature set reported by Whitfield et al that includes genes involved in cell proliferation was also examined [17]. Finally, we also assessed an ATP-binding cassette transporter (ABC) gene set that included genes involved in drug transport and was previously shown to predict chemotherapy response in cell lines by Szakacs et al [18].

Gene annotations were based on UniGene Build 185 that was used to match the genes in each of the above publication to probe sets on Affymetrix U133 A gene Chips. The gene sets are listed in Table 2 and a complete list of all probe sets that comprise each of the 852 sets is provided in supplementary Table 1 (Additional file 1).

Table 2 Gene sets used in this analysis

Gene set enrichment score was calculated as reported previously [5]. This score is a measure of the degree to which a gene set is overrepresented at the extremes of the entire ranked gene list. Significance was assessed by permuting class labels (i.e. response category) and calculating enrichment scores for the permuted data sets that yielded a null distribution. Nominal p-value for a score was derived from comparison with this null distribution. In order to adjust for multiple hypothesis testing, the false discovery rate (FDR) q-value was calculated for each gene set. The q-value could be considered as an FDR-adjusted p-value; however unlike p-values which express the probability of a false-positive result for a single test, the q-value gives an estimate of the proportion of false positives for a set of results [19]. Gene sets with FDR q-value ≤ 0.25 were considered to be of interest which indicates that the result is likely to be valid 3 out of 4 times and represents a previously proposed cut off in the literature [5]. Gene set enrichment analysis was performed using R package of GSEA (version 1.0) provided by the Broad Institute of Massachusetts Institute of Technology (Cambridge, MA).

Results

Differentially expressed genes

Twenty-seven out of 51 patients (53%) had pCR among the ER-negative tumors and 7 out of 82 (9%) among the ER-positive tumors. The much lower response rate in ER-positive cancers is consistent with previous reports in the clinical literature [1–4]. First, we examined if we could identify differentially expressed genes between cases with pCR and RD using unequal variance t-test. We performed this analysis separately for ER-negative and ER-positive cancers. Figure 1 shows results of the BUM analysis of the p values from the t-test. In ER-negative cancers, the FDR rate associated with the lowest p value (p ≤ 0.00087) was 40%. In ER-positive cancers, the FDR was close to 100% for all observed p values. The scarcity of low p-values in the ER-positive group is due to the unbalanced sample size (i.e. few informative cases, 7 pCR only) and suggests and underpowered analysis that violates t-test assumptions.

Figure 1

Distribution of p values computed from unequal variance t-test in patients with ER-negative and ER-positive tumors, respectively. A, Gene expressions were compared between ER-negative tumors that had pathologic complete response and those that had a lesser response to preoperative chemotherapy. The resulting p-values for all comparisons were modeled as beta uniform mixture. The straight line indicates the contribution of the uniform component, and the curved line is the fitted beta-distribution from the observed values. Deviation above the straight line indicate p-values that may represent true discovery. B, Distribution of p values in patients with ER-positive tumors.

These observations indicate that in these 2 data sets, t-test cannot reliably identify differentially expressed genes. However, these results do not necessarily indicate that there are no real transcriptional differences between cases with pCR compared to RD when ER-positive and ER-negative cancers are analyzed separately. It is possible, that no individual gene meets the threshold for statistical significance after correcting for multiple hypotheses testing because the transcriptional differences are modest relative to the technical noise and biological variability that is present in the data. Analysis at the single gene level may also miss small but coordinated expression differences in a larger number of genes that could belong to important biological pathways. In some situations, small coordinated change in the expression of many genes that belong to a particular metabolic pathway can have robust functional consequences [20]. Such subtle gene expression differences would not be identified easily by pair-wise comparisons using t-statistics. Different analytical tools, including gene set enrichment analysis (GSEA) were developed to test for potentially relevant but small scale transcriptional differences in pre-defined sets of genes.

Gene set enrichment analysis to identify pathways associated with complete response to preoperative chemotherapy

We applied GSEA to the 51 ER-negative tumors. Only 3 gene sets out of the 853 were enriched with an FDR q ≤ 0.25. These gene sets included; (i) the proliferation set (FDR q = 0.1, p = 0.05) (Figure 2A and Supplementary Figure 1A in Additional file 2); (ii) the genomic grade index set (FDR q = 0.04, p = 0.08) (Figure 2B and Supplementary Figure 1B in Additional file 2); and (iii) the E2F3 pathway gene set (FDR q = 0.2, p = 0.07) (Figure 2C and Supplementary Figure 1C in Additional file 2). All of these were enriched in the group with pathologic complete response, while no gene set was enriched in the group with residual cancer.

Figure 2

Gene set enrichment results for ER-negative breast cancers. Running enrichment scores (RES) and the location of each probe set within the complete ranked ordered gene list for each gene set. The first dotted line indicates the position of the maximum RES, the second dotted line indicates the zero position of the ranking metric score. A, proliferation set (probe set n = 74). B, genomic grade index (probe set n = 242). C, E2F3 pathway (probe set n = 173) (pCR = pathologic complete response; RD = residual disease). Heat maps corresponding to these plots are provided on supplementary figure 1 in Additional file 2.

We performed the same analysis on the 82 ER-positive cases. Two gene sets were enriched in the group with pathologic complete response including the proliferation set (FDR q = 0.001, p = 0.002) (Figure 3A and Supplementary Figure 2A in Additional file 3) and the genomic grade index set (FDR q = 0.015, p = 0.01) (Figure 3B and Supplementary Figure 2B in Additional file 3). In the group with residual cancer, 2 other gene sets showed enrichment including the ER-associated gene list (FDR q = 0.03, p = 0.04) (Figure 3C and Supplementary Figure 2C in Additional file 3) and the mutant p53 gene signature (FDR q = 0.0019, p = 0.07) (Figure 3D and Supplementary Figure 2 D in Additional file 3). The complete list of probes and genes included in the five enriched gene sets is presented in supplementary Table 2 (Additional file 4).

Figure 3

Gene set enrichment results for ER-positive breast cancers. Results are presented as in figure 2. A, proliferation set. B, genomic grade index: C, ER-associated genes (probe set n = 201). D, mutant p53 gene signature (probe set n = 25). Heat maps corresponding to these plots are provided on supplementary figure 2 in Additional file 3.

These results indicate that higher expression of proliferation-related genes characterized cancers with pCR among both ER-negative and ER-positive cancers. There are several single gene markers of proliferative activity; the one that is most commonly used in the clinic is Ki67 (which was also included in the proliferation gene set). We therefore examined if measuring Ki67 (MIK67) mRNA expression alone is sufficient to separate cases with pCR from those with residual disease after chemotherapy. Ki67 expression is measured by 2 distinct Affymetrix probe sets "212021_s_at" and "212023_s_at". In ER-negative cancers, neither of these probe sets were significantly differentially expressed according to pathologic response to chemotherapy (unequal variance t-test, p = 0.97 and 0.92, respectively). In ER-positive cancers, one of the 2 probe sets ("212023_s_at") showed borderline significant over-expression in the pathologic complete response group, p = 0.06.

Also, since the ER-associated gene set was enriched in ER-positive breast cancers with residual disease after chemotherapy, we examined if quantitative assessment of the ER (ESR1) mRNA alone could provide the same information [21]. There was no statistically significant difference in ER mRNA expression levels (probe set "205225_at") between cases with complete response and those with residual cancer after chemotherapy (unequal variance t-test, p = 0.09) among the ER-positive tumors.

Discussion

In the present study, we examined if we could find individual genes or gene sets that are significantly associated with extreme chemotherapy sensitivity in ER-positive and ER-negative breast cancers, respectively. These two major types of breast cancers differ in the expression of thousands of genes [6, 7, 22]. They also have substantially different sensitivity to cytotoxic treatment; ER-positive cancers are generally less sensitive to chemotherapy than ER-negative tumors [2]. Therefore, when these cancers are analyzed together, sensitivity markers tend to identify ER-negative tumors and are often dominated by genes that reflect the ER-status of the tumor [4]. In order to identify markers of response that are independent of ER-status we analyzed these 2 groups of breast cancers separately. To our surprise, the commonly used approach which performs gene wise comparison between responders and non-responders failed to identify any genes that could be declared differentially expressed with statistical confidence. The estimated false discovery rate was > 40% among the top differentially expressed genes in ER-negative cancers and the FDR was even higher among ER-positive cancers. These findings are in contrast with the results that can be obtained when the entire patient cohort is analyzed together. When we searched for differentially expressed genes including both ER-negative and -positive cases, we could identify over 400 genes with an FDR ≤ 1% [4].

We next examined if coordinated but relatively small scale differences in the expression of sets of genes that belong to functional pathways are associated with response. Such small scale differences at the individual gene level may not be readily identified by t-statistics but gene set enrichment analysis may be able to detect these. Two gene sets emerged as strongly enriched in cancers with pCR to chemotherapy in both ER-positive and ER-negative cancers. These included 44 genes (corresponding to 74 probe sets) involved in cell proliferation and 183 genes (corresponding to 242 probe sets) that distinguish histologically high grade cancers from low grade tumors. These observations are consistent with the literature that suggests that highly proliferative cancers are more sensitive to cytotoxic treatment in general [23]. However, there is no consensus how to best measure proliferative activity [24]. To underscore the power of gene set analysis, we noted that the proliferation signature as a whole was significantly overrepresented in highly chemotherapy sensitive tumors. However, a commonly used proliferation marker Ki67 that was included in the signature showed no significant over-expression when tested alone. It is also well documented in the clinical literature that high histological grade is associated with better response to preoperative chemotherapy [25]. It was reassuring to observe that the same association holds up for the genomic grade index too.

We also made 3 novel observations. Our results indicate that the expression of genes involved in the E2F3 pathway may be associated with high degree of chemotherapy sensitivity in ER-negative cancers. Considering that the E2F3 family of transcription factors plays a critical role in regulating cell cycle progression this association is not surprising [26]. Nevertheless, no previous reports linked E2F3 activity to chemotherapy response. It is also intriguing that no association between E2F3 pathway and pCR was seen in ER-positive cancers. We also observed higher expression of mutant-p53 associated genes in relatively chemotherapy resistant ER-positive breast cancers. A similar association was not seen among ER-negative cancers, which suggests that p53 dysregulation may have different consequences on chemotherapy sensitivity depending on the hormone receptor status of the cancer. This may partly explain the conflicting results about the role of p53 mutation in chemotherapy response in the literature. Some studies suggested that functional p53 defects predict for increased sensitivity to anthracycline chemotherapy [27]. Others reported that p53 mutations are associated with resistance to anthracyclines [28, 29]. It has also been shown that breast cancer cell lines exhibit different transcriptional response to chemotherapy in vitro depending on their hormone receptor status and molecular class. For example, the expression of p21, a p53-regulated protein, was highly induced in ER-positive cells, but only weakly in ER-negative breast cancer cell lines in response to anthracycline exposure [30]. This suggests that p53-mediated apoptosis may be more important in ER-positive (luminal) than in ER-negative (basal-like) cells. Future biomarker studies will need to consider the possibility that the predictive value of a biomarker may depend on the molecular subtype of the cancer [31].

We also observed that those ER-positive cancers that had low expression of ER-associated genes were more sensitive to chemotherapy. This was independent of the actual level of ER expression and indicates that some ER-positive breast cancers do not posses the full transcriptional signature of ER-activity. These cancers showed increased chemotherapy sensitivity.

Our study has limitations. All patients received combination chemotherapy that represents the current standard of care for this patient population. This makes our observations more relevant for clinical practice but at the same time limits our ability to decipher drug-specific response pathways. This could have biased our results towards detecting "generic" drug sensitivity pathways such as proliferation. The gene sets that we tested were assembled from the published literature and often contain overlapping genes represented in multiple gene sets. Our current knowledge of biology is incomplete and does not allow defining precisely all the genes that contribute to a given biological process or represent a unique molecular pathway. This may explain why the majority of the 852 gene sets that we examined including numerous apoptosis and signaling pathways did not show enrichment by chemotherapy response. It is also important to consider that GSEA is a method to demonstrate that the expression of a given gene set is overrepresented in a the top or bottom of particular gene lists ranked by correlation with clinical outcome. However, this method cannot be used to predict response in a new case. How to translate GSEA results into a prospective, single sample response predictor remains an unsolved bioinformatics challenge.

Conclusion

We found that it is difficult to identify individual genes associated with chemotherapy response with statistical confidence when ER-negative and -positive breast cancers are analyzed separately. In contrast, GSEA revealed several biological pathways that were associated with response. These included proliferation related genes and the genomic grade index that were enriched in cancers with high sensitivity to chemotherapy regardless of ER status. Genes included in the E2F3 pathway were also enriched in ER-negative and highly chemotherapy sensitive cancers. On the other hand, a mutant p53 gene expression signature and a set of highly ER-associated genes were enriched in ER-positive and chemotherapy resistant cancers. These results suggest that proliferative activity confers increased sensitivity to chemotherapy in breast cancer in general, whereas other biological pathways such p53 mutation, E2F3 activation may be more selective and influence chemotherapy sensitivity only in particular molecular subtypes of breast cancer.