A new clinical-genomic model to predict 10-year recurrence risk in primary operable breast cancer patients

This study aimed to validate the long-term prognostic value of a new clinical-genomic model, Distant Genetic Model-Clinical Variable Model 6 (DGM-CM6), developed in Asia as a prognostic panel for all subtypes of breast cancer. We included 752 operable stage I–III breast cancer patients representing all subtypes treated from 2005 to 2014 as the validation cohort. The median follow-up was 95.8 months. The low- and high-risk patients classified by DGM-CM6 (RI-DR) had significant differences in 10-year distant recurrence-free interval (DRFI) (94.1% vs. 85.0%, P < 0.0001) and relapse-free survival (RFS) (90.0% vs. 80.5%, P = 0.0003). External validation using EMTAB-365 dataset showed similar observation (P < 0.0001). DGM-CM6 was an independent prognostic factor by multivariate analysis with hazard ratios of 3.1 (1.6–6.0) for RFS (P = 0.0009) and 3.8 (1.6–9.0) for DRFI (P = 0.0028). Comparing the C-index of DGM-CM6 and PAM50-ROR scores, the former performed better than the latter in predicting long-term DRFI and RFS, especially in N0, ER/PR-positive, and HER2-negative patients.

This clinical model was based on our previous work that identified the most important prognostic factors as the number of axillary lymph nodes involved, age at diagnosis (≤40vs >40 years), prominent lymphovascular invasion (LVI), oestrogen receptor (ER) status, tumour grade, and tumour size (>2 cm) 11,12 . Incorporating both the genomic and clinical data, DGM-CM6 (recurrence index for distant metastasis [RI-DR]) proved to be the most predictive 13 .
In this study, we assessed and validated the prognostic value of DGM-CM6 (RI-DR) in different molecular subtypes of EBC after surgery based on the independent dataset.
Subgroup analyses revealed that DGM-CM6 (RI-DR) and DGM could distinguish the low-and high-risk patients in luminal, HER2, and triple-negative EBC (Supplementary Table S6). However, DGM score and RI-DR were not significant factors in patients with HER2-overexpressed and triple-negative breast cancer; the low-risk group had a trend towards a better outcome than the high-risk group. When we confined the analysis to luminal N0-N1 patients, DGM and RI-DR could significantly distinguish the low-and high-risk patients (Supplement Table S6).
For the interaction between DGM-CM6 (RI-DR) and chemotherapy, RI-DR was capable of classifying lowand high-risk N0-2 patients as 10-year DRFI regardless of chemotherapy administration. The 10-year DRFI for low-and high-risk patients who did not receive chemotherapy was 97.0% and 82.3% (P = 0.012), respectively. The corresponding rates in patients receiving chemotherapy were 93.4% and 85.2% (P = 0.0008), respectively (Fig. 1C). The int eraction between RI-DR and chemotherapy using RFI, DRFS, and RFS as study endpoints was shown in Supplementary Figs. S1-S3.
Comparison to PAM50 intrinsic subtypes. According to research-based PAM50 intrinsic subtypes, a heatmap was generated by unsupervised clustering of all 752 patients combining our genomic panel with IHC4 genes (ER, PR, HER2 and MKI67). Our gene panel differentiates each subtype correctly ( Fig. 2A); for example, BUB1B, TPX2, BLM, and DDX39 were clustered together with MKI67; furthermore, the panel was capable of distinguishing luminal A from luminal B subtypes. TRPV6 and CLCA2 were clustered together with ERBB2 and the HER2 subtype was differentiated from other subtypes.
The score distributions of DGM were significantly different among PAM50 intrinsic subtypes (Fig. 2B). Luminal A patients had the lowest DGM scores among all subtypes (P < 0.0001). PAM50 ROR score (ROR-S) low-risk patients were usually classified as low-risk by DGM; however, DGM further classified some patients in normal-like and luminal B subtypes as low-risk. DGM divided the ROR-S low-risk group into low-risk and high-risk groups; the former had a 10-year DRFI of 93.5% (86.3%, 97.0%) and the latter 77.1% (53.1%, 89.9%) (P = 0.0019) ( Table 2). DGM also identified low-risk patients in the ROR-S intermediate-risk group with a 10-year DRFI of 95.7% (87.1%, 98.6%). The gene expression levels of DGM-low and -high patients were significantly different (Supplementary Table S4).
Similarly, DGM combined with clinical variables (DGM-CM6 or RI-DR) separated ROR-S low-and intermediate-risk patients into low-and high-risk groups significantly (Table 2).
Validation in an external dataset. The performance of DGM (clinical data was inadequate to test DGM-CM6) in predicting the outcomes of N0-2 patients from the EMTAB-365 dataset revealed that the 10-year DRFS was 62.1% in the high-risk group and 82.3% in the low-risk group (P < 0.0001) (Fig. 3). According to the PAM50, the ROR-S low-, intermediate-and high-risk patients had 10-year DRFS rates of 80.1%, 67.2% and 57.8%, respectively (Fig. 3).

Discussion
The new clinical-genomic model DGM-CM6 serves as an independent prognostic factor in patients with N0-N2 primary operable breast cancer, especially the luminal subtype; however, its prognostic value in non-luminal subtype needs to be confirmed with more data. The hazard ratios for DRFI and RFS were 3.8 (1.6-9.0, P = 0.0028) and 3.1 (1.6-6.0, P = 0.0009), respectively (Table 3). This model also divided PAM50 ROR low-and intermediate-risk patients into different risk groups ( Table 2). The 10-year rates of DRFI in RI-DR low-risk and ROR low/intermediate-risk groups were excellent, ranging from 94.6% to 98.5%. The data obtained in our study suggest that our model can identify high-risk patients from the ROR low-risk group and low-risk patients from the ROR intermediate-risk group ( Table 2). As a result, 44/192 (22.9%) PAM50 luminal A patients were identified as high-risk and 43/212 (20.3%) luminal B patients as low-risk ( Table 1).
Although the multi-gene panel was initially developed without considering breast cancer subtypes, the heatmap and correlation analyses revealed that our panel can differentiate among PAM50 intrinsic subtypes ( Fig. 2A,B). The heatmap showed that the gene expression levels of BUB1B, TPX2, BLM and DDX39 are different between PAM50 luminal A and B subtypes. Other researchers have made similar observations; BUB1B is associated with poor prognosis in luminal A breast cancer 14 . TPX2 is the most well-connected gene within a proliferation network; its knockdown significantly affects metastasis but not tumour proliferation in oestrogen receptor-positive tumours 15 . Bloom syndrome helicase (BLM) has key roles in homologous recombination repair; PAM50 luminal A subtype is more likely to express low levels of BLM mRNA 16 .
Concordant statistics using the validation dataset revealed DGM-CM6 had higher C-indices than DGM and PAM50 ROR scores (Fig. 2C). This is understandable as DGM-CM6 incorporates clinical information in the model that might increase the C-index. Confined to node-negative, ER+/PR+ and HER2-negative patients, the C-indices of DGM and DGM-CM6 for DRFS and RFS were 0.72-0.75; however, the C-index of ROR-S was 0.65-0.66 (Fig. 2D). This may be related to the fact that our dataset is based on an Asian population with reduced odds of the basal-like subtype and apparent ethnicity differences 17 . The C-index of ROR-S for post-menopausal node-negative luminal women in anastrozole or tamoxifen alone or combined randomised clinical trials was reported as 0.78 18 .
The main goal of adjuvant chemotherapy is to reduce the risk of distant recurrence. The current study demonstrated very low-risk DR within 5 years in the DGM-CM6 low-risk group. However, some late recurrences developed after 5 years (Fig. 1A). Patients in the current study received hormonal therapy for only 5 years; the DR after 5 years was probably related to the duration of hormonal therapy. The type and risk of recurrence vary significantly among different molecular subtypes; furthermore, our genomic information is highly correlated with the PAM50 subtype. Numerous multi-gene panels or clinical-genomic models have been developed to assist in decision making for adjuvant systemic therapy. However, most of them focus on luminal subtypes and are rarely shown to play a role in basal-like or HER2 positive subtypes. In our gene panel, TRPV6 and CLCA2 were clustered together with ERBB2 and could differentiate HER2 from non-HER2 subtypes (Fig. 2A). Both genes are related to ion channel pathway control 19,20 TRPV6 expression leads to reduction in basal calcium influx and cellular proliferation and is significantly elevated in basal-like and HER2 subtypes 19 . CACL2 is a tumour suppressor, involved in the p53 tumour suppressor network and has a significant effect on cell migration and invasion 21 . These 2 genes could be novel targets for HR-negative breast cancer 19 .
For patients with HER2 positive breast cancer treated with curative surgery, adjuvant trastuzumab for one year is the standard care. However, identifying patients, who are at a higher risk of recurrence and would, therefore, benefit more from novel anti-HER2 agents such as pertuzumab and neratinib is paramount. There is an urgent need for a predictive tool to guide the systemic treatment strategies of these patients. Our clinical-genomic model can classify breast cancer patients into high recurrence risk and low recurrence risk regardless of molecular subtypes, which has the potential to help clinicians make more informed decisions about systemic treatments.
low-and high-risk groups 22 . Comparison of GenesWell BCT score with ODX RS revealed that BCT score classified more low-risk patients than RS in patients aged 50 years or less (73.0% versus 33.6%) 7 . Since Asian breast cancer patients are usually pre-menopausal 23 , further studies, including our model, are necessary to identify which test is more accurate in this subpopulation.
There were some limitations to our study. First, the ideal prognostic validation dataset should recruit only patients who have not received systemic therapies because the risk of recurrence after adjuvant therapy may be underestimated. We had 82 (10.9%) patients who did not receive chemotherapy, but this number was too small for further analysis. It is clear that this study cannot provide adequate information for patients to make a decision about adjuvant chemotherapy. However, the potential prognostic value of our DGM-CM6 model should be noted for the significant difference between the low-and high-risk breast cancer recurrence in large cohorts. Second, most patients with HER2-positive breast cancer did not receive anti-HER2 therapy. The utility of this model in the era of anti-HER2 treatment is unclear. Finally, only a few triple negative breast cancer patients were low-risk according to our model; further investigation is necessary for this group.
In conclusion, we developed a model combining genomic and clinical information as a prognostic tool for non-metastatic breast cancer. This multi-gene model can provide not only clinical outcome information before treatment but also may play a tool to assist in the risk-benefit judgement of systematic adjuvant treatments, especially in Asian patients.

Materials and Methods
Patient population. Breast cancer patients, who had undergone microarray analysis of their primary tumour were enrolled in this study. The Consolidated Standards of Reporting Trials (CONSORT) flow diagram for this study is shown in Fig. 4. Details of the training and testing information for DGM-CM6 has been reported in our previous publication (Supplementary Tables S1 and S2) 13 . This study focused on validation using a dataset obtained from the Affymetrix platform. www.nature.com/scientificreports www.nature.com/scientificreports/ The internal validation cohort consisted of 752 patients, who had a microarray study performed for their primary tumours. This study was performed in a prospective way that all alive participants gave written informed consent to use their frozen tumor tissues from the biobank for the purpose to identify poor or good gene expression profiling. The inclusion criteria were pathology stage pN0-2 (0-9 axillary lymph nodes were positive) breast cancer patients after primary surgery with either mastectomy or breast-conserving surgery (BCS). Patients who had preoperative chemotherapy and pN3, T4, and/or M1 disease were excluded. The protocol and informed consent documents were reviewed and approved by the institutional review board (IRB) of the Koo Foundation Sun Yat-Sen Cancer Center in Taipei, Taiwan (IRB no. 20131001 A).
The EMTAB-365 dataset was used as the external validation cohort, which is the most extensive dataset using Affymetrix U133 Plus 2.0 microarray to analyse gene expression profiles of primary tumour tissues 24 . A total of 426 patients with pN0-N2 regardless of breast subtypes and microarray data were included (http://www.ebi.ac.uk/ arrayexpress).  www.nature.com/scientificreports www.nature.com/scientificreports/ Affymetrix microarray and PAM50 subtyping. The mRNA microarray results were reported previously 8,25 RNA was extracted from primary tumour tissue using TRIZOL reagent (Invitrogen/Thermo Fisher Scientific, Waltham, MA, USA) and purified with an RNEASY Mini Kit (Qiagen, Hilden, Germany); the purity was evaluated with an Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA). According to the Affymetrix protocol, hybridisation targets were prepared from total RNA and hybridised to U133 Plus 2.0 (U133P2) arrays (Affymetrix, Santa Clara, CA, USA). The details of the study protocol were reported previously 25 . Each patient was assigned to an intrinsic molecular subtype of breast cancer (luminal A, luminal B, HER2-enriched, basal-like and normal-like) using the research-based PAM50 subtyping 26,27 .
External validation Affymetrix U133P2 dataset was obtained from ArrayExpress (EMTAB-365). Raw CEL files were pre-processed using the robust multi-array average method in the affy package of R software 28,29 . Quantile normalisation was performed to reduce potential systematic biases. The classification of PAM50 subtypes and calculation of risk of recurrence (ROR) score were performed using genefu R package 26,30,31 . Algorithm of DGM and DGM-CM6. The algorithm for the DGM is summarised as follow: Gene N, N 17 (scores rescaled to 1 100) N The RI-DR score was calculated in 2 steps: 1) the genetic score was calculated as described above; and 2) clinical and genetic scores were integrated. The algorithm is summarised as follows: