Mutated Pathways as a Guide to Adjuvant Therapy Treatments for Breast Cancer Free

Breast cancer is a highly heterogeneous disease consisting of distinct molecular subtypes having different prognostic and therapeutic responses (1). Of interest to this article is that surgery is often followed by adjuvant therapy to diminish the chance of recurrence. Chemotherapy is used to stop the growth of cancer cells by either killing them or otherwise halting division (2); hormone therapy lowers estrogen levels or blocks its action (3). The choice between the two is usually made based on the levels of (4, 5) estrogen receptor (ER), progesterone receptor (PR), and HER2 (6). In general (4), Luminal-A is treated with hormone therapy; luminal B is often treated with chemotherapy, but occasionally with hormone therapy; and HER2 and basil-like subtypes are treated with chemotherapy.

Until now the choice of therapy has been guided by the levels of RNA transcripts of the three hormone receptor genes. Our understanding of cancer biology (7) and the state of computational science (8–15) has, however, now reached a point where a much fuller profile of the cancer cell can be used to guide the choice of therapy.

The development of cancer is a complex multistep process that involves accumulation of multiple mutations that lead to dysfunction of cell signaling pathways responsible for cell growth and cell fate (16). In particular, mutations are not uniformly distributed across different cellular functions, but tend to cluster in a relatively well defined set of physiologically relevant pathways (17). Consequently, it is now generally accepted that causal mechanisms underlying transformation generally reflect the behavior of functionally coupled sets of genes (18–20). An analysis of mutational heterogeneity in the context of cellular signaling and regulatory pathways can therefore add to our understanding of cancer progression and its modulation by different therapeutic strategies (16).

We identified sets of mutated pathways that are found in a high percentage of TCGA (21) breast cancer samples derived from patients to whom adjuvant therapy had been, or was being, administered. We refer to these pathways as collaborative because they are mutated in a high percentage of the same samples, as opposed to pathways that are altered in different but overlapping sets of samples. This analysis was performed using a recently developed method, Mutational Driver Pathway Collaboration (MUDPAC; ref. 20). Our goal was to determine whether different therapy groups—the chemotherapy group (CT) and the hormone therapy group (HT)—have distinguishing sets of mutated pathways, and if they do, to determine the composition of the two groups. The central result of this article is that we found strong correlation between altered processes (collaborative pathway groups) and therapy group. This suggests a useful addition to our knowledge of cancer subtypes and provides a sound molecular basis for subtype recommendations, with perhaps some moderate shifts in recommendations. These observations could have important implications for cancer biology and therapy, which we discuss later. A flowchart of our method can be seen in Fig. 1.

Breast cancer tissue somatic mutations (.maf file) were downloaded from TCGA on March 2013. All mutations in .maf file are sequenced and annotated before any systemic treatment and therapy (22).

Tissue sampling is carried out by the TCGA project; no further tissue sampling is performed in this study. Sample subtype information was obtained from the supplementary table of ref. 21 where PAM50 is used to stratify subtypes based on mRNA data. The therapy treatment and drug information of each patient was downloaded from breast cancer TCGA clinical data released in July 2013. Only chemotherapy and hormone therapy patients are included is this study, which are the only adjuvant treatments with a sufficient number of samples in TCGA (>50). Cooperative pathways are identified using samples from patients who have been treated by chemotherapy or hormone therapy, but not both. After filtering, 56 chemotherapy samples remained (34 basal-like, 10 HER2⁺, 9 Luminal-A, 3 Luminal-B), with 3,971 mutations in 3,186 genes, and 51 hormone therapy samples remained (1 basal-like, 0 HER2⁺, 40 Luminal-A, 10 Luminal-B), with 2,049 mutations in 2,049 genes.

A total of 269 pathways (.xml and .gene files) were downloaded from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database released on June 2013; 200 of them are used in this study after excluding 3 global metabolic pathways and 66 human disease pathways.

MUDPAC is a method that identifies mutated collaborative pathway groups, i.e., pathways that are altered in a high percentage of the same samples of a certain cancer or cancer subtype (20). It consists of two steps. The first step is an identification of pathways showing statistically significant differences between nonsynonymous mutation group and synonymous mutation group, based on some characteristics of mutations, such as mutation frequency, functional mutation score, mutual exclusivity, and network topology. The second step uses a greedy search for collaborative pathways. Two more criteria are considered when selecting a new pathway into collaborative set. First, the newly selected pathway along with all pathways already in collaborative set should have a maximal coverage rate (MCR) higher than the highest mutation rate of genes in this new pathway by a given threshold (3% in this study). Second, permutation is used to test whether MCR of this new pathway is significantly higher than a background MCR (P < 0.05 in this study).

Each of the missense mutations in .maf file is assigned a functional mutation score by MutationAssessor (23) based on hg19 reference genome. Because MutationAssessor can only evaluate missense polymorphisms, the functional scores for the remaining mutations are assessed following the same criteria in MUDPAC: the highest score that can be calculated using MutationAssessor is assigned to all indels, nonsense mutations, and splice site mutations; the lowest score from MutationAssessor is allocated to synonymous mutations; the average score of all missense mutations is assigned to other missense mutations that cannot be calculated using MutationAssessor.

MuSiC (24) is used to calculate the total number of bases for each gene having available alignment data, as part of the input required by MUDPAC. MuSiC uses Broad Institute's analysis infrastructure Firehose to count bases with sufficient coverage of each gene from the given wiggle tract format file. Wiggle files contain dense, continuous data, such as GC percent, probability scores, and transcriptome data, and were downloaded from Broad Institute Firehose webpage on August 2012. The thresholds for sufficient coverage are at least 8-fold read depth in normal tissue, and at least 14-fold read depth in cancer tissue.

The Fisher exact test is used to identify KEGG pathways that are statistically enriched in mutated genes of identified pathways compared with the human genome background, and to assess correlations between pathway groups and therapy.

Altered physiologic functions in CT and HT groups were examined separately using MUDPAC. We identified three collaborative pathways (Fig. 2A) in CT, having an MCR of approximately 79%, which means these three pathways are altered in the same 44 (of 56) chemotherapy patients. Of the 51 patients who received only hormone therapy, 20 collaborative pathways (Fig. 2B) were identified with an MCR of 57%, which results in a much longer and more stable plateau. None of these 20 are among the three pathways that form the CT group (Fig. 2C).

All three mutated pathways (PI3K-Akt signaling, p53 signaling, Wnt signaling) in the chemotherapy pathway group (CTPG) are TP53 related. This is consistent with the sample distribution within this therapy group, which has 34 (61%) basal-like samples dominated by the mutated TP53 suppressor gene (21), which is known to be associated with the basal-like subtype.

In contrast with these three pathways, which play key roles in signal transduction and cell growth (Fig. 2D), the 20 pathways in the hormone therapy pathway group (HTPG) tend to have an organismal systems classification (Fig. 2D), including the endocrine, immune, central nervous, development, and digestive systems. More than half are either immune system or endocrine system related: 30% for the former (six pathways: Fc gamma R–mediated phagocytosis, chemokine signaling, B-cell receptor signaling, Toll-like receptor signaling, T-cell receptor signaling, leukocyte transendothelial migration) and 20% for the latter (four pathways: insulin signaling, estrogen signaling, progesterone-mediated oocyte maturation, prolactin signaling).

The gene sets that trigger pathway dysfunction in the two therapy groups are largely distinct from one other. The top 10 genes with highest mutation frequency in CT and HT are shown in Table 1.

Of the 10 mutated genes in each group that occur with the highest frequency, 3—TP53, PIK3CA, OBSCN—are common to the two groups. TP53 ranks first with a frequency of 75% in the CT but 10th (frequency of 8%) in HT. PIK3CA ranks first in HT with a frequency of 53%, but has a frequency of only 18% in CT. OBSCN has comparable mutation rates in the two groups. However, the mutation type and functional mutational score for OBSCN are different in the two groups. OBSCN has 10 mutations (one silent, six missense mutations, one frame shift mutation, two nonstop mutations) in CT with an average mutation score of 1.84, but it only has five mutations (two silent mutations, two missense mutations, one nonsense mutation) in HT with an average functional score of 0.44. This suggests that comparable frequencies do not necessarily imply comparable physiologic impact.

We also compared all the genes in the identified pathways of the two therapy groups with 87 plausible breast cancer genes (25): 37 confirmed drivers, and another 50 driver candidates identified with high likelihood computationally. Of the 532 genes in the CTPG, 27 are among the 87 plausible drivers (Fisher exact test P < 10⁻¹⁵). Of the 1,059 genes in the HTPG, 40 are among the 87 plausible drivers (Fisher exact test P < 10⁻¹⁵). These results suggest that the mutated genes in the identified pathways are highly enriched in drivers. For the two sets of drivers—27 genes for CTPG and 40 for HTPG—21 are common. The remaining 6 (19) are unique to the chemo (hormone) therapy patients.

The choice of adjuvant therapy, as noted above, generally follows subtype (4): the Luminal-A subtype is treated with hormone therapy; the Luminal-B subtype with chemotherapy, but occasionally with hormone therapy, and the HER2 and basil-like subtypes with chemotherapy.

The patient population in this analysis deviates somewhat from these guidelines, showing 9 of the 49 luminal-A patients are in chemotherapy. The deviation for luminal-B is more severe, with only 3 of the 13 in chemotherapy. We will comment on this below.

The subtypes and their therapy recommendations are approximately correlated. Thus, for basal-like subtype patients having chemotherapy, 95 of 102 pathways (34 basal-like patients × 3 identified pathways in chemotherapy) are mutated in CT while 7 are not; 264 out of 680 (34 basal-like patients × 20 identified pathways in hormone therapy) pathways are mutated in HT while 416 are not. Consequently, basal-like patients having chemotherapy are significantly more likely to be mutated in CTPG than in HTPG (P < 10⁻¹⁵). Similar results are observed for HER2⁺ (P = 7.35e–4) and Luminal-A (P = 2.10e–8), which are recommended for chemotherapy and hormone therapy, respectively. But there are no significant correlations found in Luminal-B patients, no matter in CT (P = 0.6119) or in HT (P = 0.7489).

There are, as expected, strong correlations between the mutated pathway groups seen in patients and the kind of therapy they receive. The distribution of mutations between HTPG and CTPG for CT patients is highly significant. In particular, the probability that the observed distribution of mutations between the two pathway groups due to chance, given the patients are in chemotherapy (Fig. 3A), is P < 10⁻¹⁵. This follows from a Fisher exact test, based on the observation that for CT (i) 141 of 168 (56 × 3) CTPG pathways are mutated and 27 are not; and (ii) 503 of 1,120 (56 × 20) HTPG pathways are mutated and 617 are not (Fig. 3A; CT columns).

Similar results apply to patients receiving hormone therapy (Fig. 3A; HT columns): 69 of 153 (51 × 3) CTPG pathways are mutated while 84 are not, and 677 of 1,020 (51 × 20) HTPG pathways are mutated while 343 of them not. The probability that mutations are equally likely in the two pathway groups, given that patients are in hormone therapy, is P = 4 × 10⁻⁷.

These results provide a previously unnoticed molecular rationale for the choice of therapy based on cancer subtype. In addition, as discussed below, the structure of Fig. 3A suggests well-defined stratifications within the standard subtypes. The implications of these observations for therapeutic choices are discussed below.

One of the three CTPG pathways—phosphatidylinositol 3′–kinase (PI3K)-Akt—is altered in almost all patients, whereas the other two are mutated in a high percentage of basal-like (32/35) and HER-2⁺ (9/10) patients, and only occasionally in luminal-A (5/49) and luminal-B (4/13) patients, irrespective of whether they received hormone or chemotherapy therapy. Consequently, these two pathways discriminate extremely well between HER2⁺ and basal-like subtypes on the one hand, and Luminal-A and Luminal-B subtypes on the other (Fisher exact test P < 10⁻¹⁵). Because the first two subtypes received chemotherapy, and rarely hormone therapy, these pathway groups also discriminate between the two treatment groups.

The situation for the HTPG is more complicated. Basal-like and HER2⁺ are typically treated with chemotherapy, and the basal-like samples do show a relatively low frequency of mutated HTPG pathways. For basal-like patients, therefore, the frequency of HTPG pathways mutated is a good correlate of treatment of choice. We also find that 29 of the 49 Luminal-A patients have virtually all HTPG pathways mutated. Three of the 9 Luminal-A patients in chemotherapy are among them. This suggests potential utility of the pathway groups, because the standard recommendation for Luminal-A patients is hormone therapy. Of the 40 Luminal-A patients who received hormone therapy, 26 have mutated HTGP pathways, so again there is reasonable correlation between the frequency of mutated pathway groups and recommended treatment. This is seen strikingly in Fig. 3A.

For HER2⁺, however, there appears to be little if any correlation. Half the HER2⁺ patients (all of whom are in chemotherapy) have mutations in almost all pathways in HTPG, and half have mutations in only a small percentage of pathways. There also appears to be little or no relation between HTPG pathways and therapy group for Luminal-B patients, nor does the assignment of Luminal-B to therapy group conform to subtype classification (10 in HT and 3 in CT, whereas guidelines would have it roughly reversed). Consequently, the classification of Luminal-B sheds no light on HTPG.

The Luminal-A pattern provides strong statistical evidence that the population of Luminal-A patients in this study can be split into two distinct groups: 29 patients have almost all 20 pathways from HTPG mutated (577 mutated, 3 not), whereas 17 patients show relatively infrequent mutations (65 mutated, 275 not). The probability that this pattern is the result of chance is P < 10⁻¹⁵. This again suggests that traditional subtypes can be stratified further and that higher resolution can be identified by pathway groups. Similarly, although the HER2⁺ pattern does not support correlation with traditional clinical guidelines for therapy, it does suggest two distinct subpopulations.

We demonstrated that HT and CT samples are characterized by mutually exclusive sets of mutated pathways that are strong correlates of subtype-based therapy recommendation. Furthermore, genes in these altered pathways are highly enriched in breast cancer drivers and may provide further guidelines for the personalized therapy. Because all the mutations reported in this article are from DNA sequenced and annotated prior to surgery or therapy, using tissue from the primary tumor at the initial site of the cancer, the results are not a consequence of mutations that might be induced by therapy.

The discovery that the altered cellular processes that drive transformation in patients having hormone therapy have a substantially different biology than those that do not may open an opportunity for the identification of new therapies. Particularly, we found that 100% of the pathways in chemotherapy are TP53 related, whereas 70% of pathways in hormone therapy are involved in organismal system, especially in immune system (30%) and endocrine system (20%).

Both immune (26, 27) and endocrine systems (28–30) can work primarily through ERs, regulating cells in pathways from these two systems could therefore lead to the successful development of therapeutic concepts in breast cancer. In brief, we hope our identified pathways could be considered as a signature to classify patients for their treatments in the future, as an auxiliary of current protocol, after we could get enough samples to test. Moreover, speaking of each individual therapy option, like hormone therapy for example, beyond targeting estrogen-related genes to lower or stop estrogen level, our mutated pathways could also open up many other treatment options by targeting other mutations that up or downregulate estrogen-related genes or initially trigger the cancer.

The mutated genes in our identified pathways tend to be sample specific (Fig. 3B). Except for TP53 and PIK3CA, which have the highest mutated rates in chemotherapy (75% in CT and 8% in HT) and hormone therapy (53% in HT and 18% in CT) groups, respectively, almost all the other genes are dispersed among samples, even for subtypes that have similar pathway patterns (i.e., have a significantly large number of mutated pathways in common). For example, 4 of the HER2⁺ patients (TCGA-A2-A04W, TCGA-A2-A0D1, TCGA-A2-A0EQ, TCGA-AO-A0JE) and 3 of the Luminal-B patients (TCGA-A2-A0CW, TCGA-A2-A0SW, TCGA-E2-A15K) have all collaborative pathways mutated (3 in CT, 20 in HT), i.e., they have identical mutation patterns. But the mutated genes in these pathways are specific to each patient: there are totally 12 genes mutated in HER2⁺ (PIK3CA, TP53, ABL2, ATR, NKD2, TXK, XCR1, ADCY8, FLT4, PPP3CA, STAT5B, IRF9), total 15 genes mutated in Luminal-B (PIK3CA, TP53, RPS6KA6, TCF7, FIGF, PAK7, PIK3R1, BCL2L11, CCNB1, FGFR2, FN1, MYH10, NFATC4, PIK3R3, TRIP10). Besides PIK3CA, which is mutated in all 7 patients, and TP53, which is mutated in 6 of 7 patients, all the remaining genes are mutated in only 1 patient, and no two genes are mutated in the same patient.

The conventional way to determine whether to apply chemotherapy or hormone therapy is according to the breast cancer molecular subtypes, which are induced by gene regulation level of several molecular factors like ER, PR, and human epidermal growth factor. This strategy has been developed for decades and greatly helped physicians to draw up clinical schemes, which are also evidenced by our identified mutated pathways: different molecular subtypes tend to have distinct mutation patterns and are approximately concordant with pathways of their recommended therapy groups.

The subtype information, including biologic process and molecular functions in epigenetic alone, however, may not fully uncover the whole picture. Mutated drivers, including both genes and pathways that initiate cancer development, may supplement this picture. We present in Results that even for the same subtype like Luminal-A, the mutational landscape is diverse between patients in CT and HT. This may advise that current treatment protocol need to be personalized by integrating genetic factors like driver genes and driver pathways.

One limitation of this article is what we reported are exploratory results only. We discussed summarized clinical information from TCGA with observation and detailed interpretation, but without very strong statistical support due to small portion of available samples in TCGA (around 50 samples for each therapy), and without further verification in new samples due to limited existing clinical sources.

Taken together, our findings suggest that systematic mutation analysis of breast cancer can reveal pathway dysfunction status and mutated gene portraits that may not easily be discovered in transcriptome level. This brings new insights about personalized treatment and targeted agents.

No potential conflicts of interest were disclosed.

Conception and design: Y. Liu, Z. Hu

Development of methodology: Y. Liu, Z. Hu

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): Y. Liu, C. DeLisi

Writing, review, and/or revision of the manuscript: Y. Liu, Z. Hu, C. DeLisi

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): Y. Liu, Z. Hu

Study supervision: Z. Hu, C. DeLisi

This work was supported by NIH grant R01GM103502-05 to C. DeLisi.

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.