Predicting Outcomes of Hormone and Chemotherapy in the Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) Study by Biochemically-inspired Machine Learning

Genomic aberrations and gene expression-defined subtypes in the large METABRIC patient cohort have been used to stratify and predict survival. The present study used normalized gene expression signatures of paclitaxel drug response to predict outcome for different survival times in METABRIC patients receiving hormone (HT) and, in some cases, chemotherapy (CT) agents. This machine learning method, which distinguishes sensitivity vs. resistance in breast cancer cell lines and validates predictions in patients; was also used to derive gene signatures of other HT (tamoxifen) and CT agents (methotrexate, epirubicin, doxorubicin, and 5-fluorouracil) used in METABRIC. Paclitaxel gene signatures exhibited the best performance, however the other agents also predicted survival with acceptable accuracies. A support vector machine (SVM) model of paclitaxel response containing genes ABCB1, ABCB11, ABCC1, ABCC10, BAD, BBC3, BCL2, BCL2L1, BMF, CYP2C8, CYP3A4, MAP2, MAP4, MAPT, NR1I2, SLCO1B3, TUBB1, TUBB4A, and TUBB4B was 78.6% accurate in predicting survival of 84 patients treated with both HT and CT (median survival ≥ 4.4 yr). Accuracy was lower (73.4%) in 304 untreated patients. The performance of other machine learning approaches was also evaluated at different survival thresholds. Minimum redundancy maximum relevance feature selection of a paclitaxel-based SVM classifier based on expression of genes BCL2L1, BBC3, FGF2, FN1, and TWIST1 was 81.1% accurate in 53 CT patients. In addition, a random forest (RF) classifier using a gene signature ( ABCB1, ABCB11, ABCC1, ABCC10, BAD, BBC3, BCL2, BCL2L1, BMF, CYP2C8, CYP3A4, MAP2, MAP4, MAPT, NR1I2,SLCO1B3, TUBB1, TUBB4A, and TUBB4B) predicted >3-year survival with 85.5% accuracy in 420 HT patients. A similar RF gene signature showed 82.7% accuracy in 504 patients treated with CT and/or HT. These results suggest that tumor gene expression signatures refined by machine learning techniques can be useful for predicting survival after drug therapies.


Amendments from Version 2 Introduction
Current pharmacogenetic analysis of chemotherapy makes qualitative decisions about drug efficacy in patients (determination of good, intermediate or poor metabolizer phenotypes) based on variants present in genes involved in the transport, biotransformation, or disposition of a drug. We have applied a supervised machine learning (ML) approach to derive accurate gene signatures, based on the biochemically-guided response to chemotherapies with breast cancer cell lines 1 , which show variable responses to growth inhibition by paclitaxel and gemcitabine therapies 2,3 . We analyzed stable 4 and linked unstable genes in pathways that determine their disposition. This involved investigating the correspondence between 50% growth inhibitory concentrations (GI 50 ) of paclitaxel and gemcitabine and gene copy number, mutation, and expression first in breast cancer cell lines and then in patients 1 . Genes encoding direct targets of these drugs, metabolizing enzymes, transporters, and those previously associated with chemo-resistance to paclitaxel (n=31 genes) were then pruned by multiple factor analysis (MFA), which indicated that expression levels of genes ABCC10, BCL2, BCL2L1, BIRC5, BMF, FGF2, FN1, MAP4, MAPT, NKFB2, SLCO1B3, TLR6, TMEM243, TWIST1, and CSAG2 could predict sensitivity in breast cancer cell lines with 84% accuracy. The cell line-based paclitaxel-gene signature predicted sensitivity in 84% of patients with no or minimal residual disease (n=56; data from 5). The present study derives related gene signatures with ML approaches that predict outcome of hormone-and chemotherapies in the large METABRIC breast cancer cohort 6 .

SVM (Support Vector Machine) learning:
Previously, paclitaxelrelated response genes were identified from peer-reviewed literature, and their expression and copy number in breast cancer cell lines were analyzed by multiple factor analysis of GI 50 values of these lines 2 (Figure 1). Given the expression levels of each gene, a SVM is evaluated on patients by classifying those with shorter survival time as resistant and longer survival as sensitive to hormone and/or chemotherapy using paclitaxel, tamoxifen, methotrexate, 5-fluorouracil, epirubicin, and doxorubicin. The SVM was trained using the function fitcsvm in MATLAB R2014a 7 and tested with either leave-one-out or 9 fold cross-validation (indicated in Table 1). The Gaussian kernel was used for this study, unlike Dorman et al. 1 which used the linear kernel. The SVM requires selection of two different parameters, C (misclassification cost) and sigma (which controls the flexibility and smoothness of Gaussians) 8 ; these parameters determine how strictly the SVM learns the training set, and hence if not selected properly, can lead to overfitting. A grid search evaluates a wide range of combinations of these values by parallelization. A Gaussian kernel selects the C and sigma combination that lead to the lowest cross-validation misclassification rate. A backwards feature selection (greedy) algorithm was designed and implemented in MATLAB in which one gene of the set is left out in a reduced gene set and the classification is then assessed; genes that maintain or lower the misclassification rate are kept in the signature. The procedure is repeated until the subset with the lowest misclassification rate is selected as the optimal subset of genes. These SVMs were then assessed for their ability to predict patient outcomes based on available metadata (see Figure 1 and reference 1). Interactive prediction using normalized expression values as input is available at http://chemotherapy.cytognomix.com.
RF (Random Forest) learning: RF was trained using the WEKA 3.7 9 data mining tool. This classifier uses multiple random trees for classification, which are combined via a voting scheme to make a decision on the given input gene set. A grid search was used to optimize the maximum number of randomly selected genes for each tree in RF, where k (maximum number of selected genes for each tree) was set from 1 to 19. Figure 2 depicts the therapy outcome prediction process of a given patient using a RF consisting of a series of decision trees derived from different subsets of paclitaxel-related genes.
Augmented Gene Selection: The most relevant genes (features) for therapy outcome prediction were found using the Minimum Redundancy and Maximum Relevance (mRMR) approach 10 . mRMR is a wrapper approach that incrementally selects genes by maximizing the average mutual information between gene expression features and classes, while minimizing their redundancies: where f i corresponds to a feature in gene set S, I(f i ,C) is the mutual information between f i and class C, and I(f i ,f j ) is the mutual information between features f i and f j .
Performance was evaluated with several metrics. WEKA determined accuracy (ACC), the weighted average of precision and F-measure, the Matthews Correlation Coefficient (MCC) and the area under ROC curve (AUC). The predicted and expected response to treatment for each individual METABRIC patient for each analyses listed in Table 1,  Table 2 and Table 3 are indexed. Patients sensitive to treatment are labeled with '0' while resistant patients are labeled '1'.

Results and discussion
The performances of several ML techniques have been compared such that they distinguish paclitaxel sensitivity and resistance in METABRIC patients using its tumour gene expression datasets. We used mRMR to generate gene signatures and determine which genes are important for treatment response in METABRIC patients. The paclitaxel models are more accurate for prediction of outcomes in patients receiving HT and/or CT compared to other patient groups.
SVMs and RF were trained using expression of genes associated with paclitaxel response, mechanism of action and stable genes in the biological pathways of these targets (Figure 3). Pair-wise comparisons of these genes with those from MammaPrint and Oncotype Dx (other genomic classifiers for breast cancer) find that these signatures are nearly independent of each other, with only a single gene overlap. The distinct differences of these signatures are due to their methodology of derivation, based on different principles and for different purposes (i.e. drug response for a specific reagent). SVM models for drugs used to treat these patients were derived by backwards feature selection on patient subsets stratified by treatment or outcome (Table 1). The highest SVM accuracy was found for the paclitaxel signature in patients treated with HT and/or adjuvant chemotherapy (78.6%). Since some CT patients were also treated with tamoxifen, methotraxate, epirubicin, doxorubicin and 5-fluorouracil, we also evaluated the performance of models developed for these drugs using the same algorithm. These gene signatures also had acceptable performance (accuracies between 71-76%; AUCs between 0.686 -0.766). Leave-one-out validation (CT and HT, no treatment, and deceased patients) exhibited higher model performance than 9-fold crossvalidation (CT and/or HT, including patients treated with radiation).
The RF classifier was used to predict paclitaxel therapy outcome for patients that underwent CT and/or HT ( Table 2). The best performance achieved with RF showed an 85.5% overall accuracy using a 3-year survival threshold for distinguishing therapeutic resistance vs. sensitivity for those patients that underwent HT.
The best overall accuracy and AUC (sensitivity and specificity) for CT/HT patients using mRMR feature selection for SVM predicting outcome of paclitaxel therapy was obtained for CT patients with 4-year survival (Table 3). Outcomes for HT patients with   . RF decision tree diagram depicts the therapy outcome prediction process of a given patient, using a RF consisting of k decision trees. Several DTs are built using different subsets of paclitaxel-related genes. The process starts from the root of each tree and if the expression of the gene corresponding to that node is greater than a specific value, the process continues through the right branch, otherwise it continues through the left branch until it reaches a leaf node; that leaf represents the prediction of the tree for that specific input. The decisions of all trees are considered and the one with the largest number of votes is selected as the patient outcome.  3-year survival were predicted with 85.7% accuracy; however, the specificity was lower in this group. SVM combined with mRMR further improved accuracy of feature selection and prediction of response to hormone and/or chemotherapy based on survival time than either SVM or RF alone. Predicted treatment responses for individual METABRIC patients using the described ML techniques are indicated in Dataset 1.
To assess model co-variate accuracy, predictions described in Table 1-Table 3 were broken down by subtype (available in Supplementary file 1). Subtypes with <20 individuals for a particular treatment combination were not analyzed. The deviation in classification accuracy between subtypes was mostly consistent with the average. One exception involved the RF and mRMR analyses, which was 8.3 to 23.0% below the average for (ER)negative, (HER2)-positive and basal subtypes in patients treated with HT. However, this deviation was not observed for CT-treated patients with the (ER)-negative subtype, which was consistent with the fact that CT response was derived from the paclitaxel gene set. (ER)-negative patients primarily received CT 6 . Further, the accuracy of the SVM models tested with CT and HT-treated patients was significantly higher for (HER2)-positive patients (26 correct, 3 misclassified; 90% accurate) compared to (HER2)negative patients (40 correct, 15 misclassified; 73% accurate).
MAPT expression (present in reduced 'CT and HT' paclitaxel For patients treated with CT with ≥4 Yr survival and CT+ HT for ≥ 5 Yr, the cost for the mRMR model was set to 64. Of those treated with CT for ≥ 4 Yr, genes were selected using a greedy, stepwise forward search, while in other cases, greedy stepwise backward search was used. Also, gamma = 0 in all cases. We assessed the separate Discovery and Validation datasets, respectively, as training and test sets and repeated the previous experiments. In this scenario, the performance of the model was poor (slightly better than random). This occurred because the gene expression distributions of many of the paclitaxel-related genes in our signature were not reproducible between these two sets (based on Wilcoxon rank sum test, Kruskal-Wallis test and t-tests; Supplementary file 2). Cross-study validation allows for the comparison of classification accuracy between the generated gene signatures. The observed heterogeneity in gene expression highlights one of the many challenges of cross-validation of gene signatures between these data from the same study exhibit drastic differences (for example, BCL2L1; Supplementary file 2). Furthermore, these gene expression differences also affect the performance of these methods when these datasets were combined (compare Table 2 and Table 4 for RF; Table 3 and Table 5 for mRMR). We considered the possibility that the Discovery model might be subject to overfitting. We therefore performed cross-study validation of the Discovery set-signature with an independentlyderived dataset (319 invasive breast cancer patients treated with paclitaxel and anthracycline chemotherapy 5 ). The mRMR+SVM CT-models performed well (4-year threshold model had an overall accuracy of 68.7%; 3-year threshold model exhibited lower overall accuracy [52%], but was significantly better at predicting patients in remission [74.2%]).
To evaluate the paclitaxel models without relying on the Validation dataset, the Discovery set was split into two distinct parts, consisting of 70% of the patient samples randomly selected for training, and a different set of 30% of samples for testing. This procedure was repeated 100 times using different combinations of training and test samples, and the median performance of these runs is reported (Table 4 and Table 5). We also compared the performance of our mRMR+SVM model with the K-TSP model 12 ( In addition, we evaluated the performance of the model derived from the Discovery set on a different set of patients treated with paclitaxel 5 . These results suggest that the aforementioned issue with Discovery training and Validation testing was primarily due to a batch effect, rather than to overfitting. While not a replication study sensu stricto, the initial paclitaxel gene set used for feature selection was the same as in our previous study 1 . Predictions for the METABRIC patient cohort, which was independent of the previous validation set 5 used in Dorman et al. 1 , of the either same (SVM) or different ML methods (RF and SVM with mRMR) exhibited comparable or better accuracies than our previous gene signature 1 .
These techniques are powerful tools which can be used to identify genes that may be involved in drug resistance, as well as predict patient survival after treatment. Future efforts to expand these models to other drugs may assist in suggesting preferred treatments in specific patients, with the potential impact of improving efficacy and reducing duration of therapy.

Conclusion
In this study we used METABRIC dataset to predict outcome for different survival times in patients receiving hormone (HT) and, in some cases, chemotherapy (CT) agents. We used published literature and various machine learning methods in order to identify optimal subsets of genes from a biologically-relevant initial gene set that can accurately predict therapeutic response of patients who have received chemotherapy, hormone therapy or a combination of both treatments. The SVM methodology has been previously shown to outperform randomized gene sets 1 . The predictions made by our method are based on the level of an individual drug. Genomic information has been shown to correlate with tumor therapy response in previous studies 5,13-17 . From these studies, analytical methods have been used to develop gene signatures for chemotherapy resistance prediction 5 , subtypes (PAM50), and metastatic risk stratification (Oncotype DX™, MammaPrint ® ). We also examined the method exhibiting the best performance in the Sage Bionetworks / DREAM Breast Cancer Prognosis Challenge 18 , which was also phenotype-based, however it produces outcome signatures based on molecular processes rather than the cancer drugs themselves. While interesting and informative, the results cannot be directly compared. Our approach may be useful for selecting specific therapies in patients that would be expected to produce a favorable response. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Supplementary File 1: Accuracy of SVM, RF and mRMR by Patient Co-variates.
Accuracy of all models described in Table 1-Table 3 were further broken down by provided patient subtype information (ER, HER, PR, PAM50, and LN).
Click here to access the data.

Supplementary File 2: Variation of Gene Expression Distribution between Discovery and Validation Datasets.
Whisker plots showing the distribution of expression in the Discovery and Validation METABRIC datasets for 26 genes used in the paclitaxel gene signature.
Click here to access the data. 1.

1.
Thank you for your suggestion. As recommended, we have repeated the 70/30% cross-validation analysis performed in the manuscript (Tables 4 and 5) with the same genes obtained from the dataset (Tables 4 and 5), but using the dataset alone for training and testing. Discovery Validation We found that this analysis had a similar performance level as the analysis reported in the main manuscript (Tables 4 and 5). There are exceptions. The mRMR+SVM gene signature developed using "CT-only" patients at a 5-year threshold was much less accurate using the data. Validation However, the "CT-only" subset of the dataset is small (N=31), and thus variability is not Validation unexpected. Overall, this analysis suggests that the cross-validation issue was indeed mostly due to batch effects.
The following sentence was written in the main text which describes this result: "Starting with the same set of genes, we also trained a separate model using the The authors were very responsive to the previous round of reviews, including more robust cross-validation and cross-study validation and comparison with other classifiers. Particular concerns remain that the author's conclusions that it is inappropriate to perform cross-study validation due to batch effects are incorrect, particularly since this challenging task is essential to assess overfitting and for clinical translation of classifiers. In addition, the conclusion was insufficiently revised to place their classifier in the context of the broader literature in this field.

Methods
Abbreviations SVM and RF must be spelled out as Support Vector Machine and Random Forest on first use. This was not addressed in the revised methods section.

Results
The authors did perform a robust cross-study validation, as requested in the previous review. We agree this is challenging, due in part to batch effects as reported in this manuscript. However, such cross-study validation is essential to assess the accuracy of classifiers. It is also essential to have translation of genomic signatures into the clinic, where even different assays may be used. To address these concerns the authors must do the following: (a) Remove the sentence "This heterogeneity indicates that it is inappropriate to test our gene expression signatures derived by one of these datasets using the other dataset." (b) Discuss the importance of cross-study 1.

2.
heterogeneity indicates that it is inappropriate to test our gene expression signatures derived by one of these datasets using the other dataset." (b) Discuss the importance of cross-study validation, challenges in this application, and potential of overfitting of suggested by these results.
The author's response that specific therapies were not provided in METABRIC is incorrect. According to Curtis , (2012) "Nearly all oestrogen receptor (ER)-positive and/or lymph node et al. (LN)-negative patients did not receive chemotherapy, whereas ER-negative and LN-positive patients did. Additionally, none of the HER2 patients received trastuzumab. As such, the treatments were homogeneous with respect to clinically relevant groupings." Therefore, the previous criticism #12 remains. Covariates such as ER/HER2/LN or PAM50 subtypes must be included in a table describing the sample cohorts remains. In addition, accuracy must be computed separately for these co-variates or included in the machine learning model.

Conclusion
The discussion is insufficient. It still lacks sufficient context of existing genomics classifiers in the literature. The discrepancy between their algorithm and clinical assays is confusing in revised sentence "Unlike Mammaprint and Oncotype Dx tests, this model focuses on predicting survival prediction based on gene expression in the tumor, presumably before or during drug therapy." As written, it appears to disregard the long history of predicting clinical outcome from gene expression involved in developing these classifiers from gene expression data (e.g., van't Veer 2002) et al., into clinical assays based upon expression of smaller numbers of genes.
Based on the previous review, the authors include context with other predictions of the METABRIC data in the response to the reviewers. This must also be included in the Conclusion to assess the relevance of their findings in the literature. Response: These abbreviations are now spelled out upon their first use in the main text (Methods section).

Results
Comment 2. The authors did perform a robust cross-study validation, as requested in the previous review. We agree this is challenging, due in part to batch effects as reported in this manuscript. However, such cross-study validation is essential to assess the accuracy of classifiers. It is also essential to have translation of genomic signatures into the clinic, where even different assays may be used. To address these concerns the authors must do the following: (a) Remove the sentence "This heterogeneity indicates that it is inappropriate to test our gene expression signatures derived by one of these datasets using the other dataset." (b) Discuss the importance of cross-study validation, challenges in this application, and potential of overfitting of suggested by these results.
Response: In regards to this point: (a) This sentence has been removed, as requested. (b) To address concerns regarding potential overfitting of our models, we cross-validate the acquired models to a non-METABRIC data set (from an independent study). In the Sage Bionetworks / DREAM Breast Cancer Prognosis Challenge, cross-study validation was performed using the "OsloVal" data set, which consists of gene expression and copy number data from 184 breast cancer patients (Margolin , 2013). However, this dataset is not publically available and et al. requires Ethics Board / IRB Review which we did not believe to be worth the effort. Instead, we performed cross-study validation on the gene expression of 310 breast cancer patients made publically available by Hatzis (2011). et al.
Analysis of this dataset was successful for the mRMR + SVM models developed using chemotherapy-treated patient ("CT" models), where the threshold for resistance was set to 3-years and 4-years. The "CT 3-year" model performed well predicting responsive patients (74.2% accuracy), while the "CT 4-year" model performed better predicting non-responsive patients (75.1% accuracy). The "CT 4-year" model outperformed the "CT 5-year" model for both sensitive and resistant patient data sets.
Random Forest and mRMR+SVM models which used hormone-treated patients ("HT" and "CT+HT") were much less accurate compared to the "CT-only" models, and predict patients a large percentage of patients from the Hatzis data as sensitive.
In the main manuscript, we have replaced the removed sentence from (a) and have written the following: "Cross-study validation allows for the comparison of classification accuracy between the generated gene signatures. The observed heterogeneity in gene expression highlights one of the many challenges of cross-validation of gene signatures between these data from the same study exhibit drastic differences (for example, ; Supplementary file 2). Furthermore, these gene BCL2L1 expression differences also affect the performance of these methods when these datasets were combined (compare Table 2 and Table 4 for RF; Table 3 and Table 5 for mRMR). We considered the possibility that the Discovery model might be subject to overfitting. We therefore performed cross-study validation of the Discovery set-signature with an independently-derived dataset (319 invasive breast cancer patients treated with paclitaxel and anthracycline chemotherapy ). The mRMR+SVM CT-models performed well (4-year threshold model had an overall accuracy of 1.

2.
mRMR+SVM CT-models performed well (4-year threshold model had an overall accuracy of 68.7%; 3-year threshold model exhibited lower overall accuracy [52%], but was significantly better at predicting patients in remission [74.2%])." Additionally, none of the HER2 patients received trastuzumab. As such, the treatments were homogeneous with respect to clinically relevant groupings." Therefore, the previous criticism #12 remains. Covariates such as ER/HER2/LN or PAM50 subtypes must be included in a table describing the sample cohorts remains. In addition, accuracy must be computed separately for these co-variates or included in the machine learning model.

References
Response: Thank you for the clarification regarding patient treatment. As a response, we have added an additional supplementary table which breaks down the accuracy of our models by subtype (ER, HER2, PR, LN and PAM50; Dataset 2). In the main text, we note that accuracy of most models are consistent between subtypes (+/-10% deviation in accuracy). Subtypes with less than twenty individuals were ignored due to its small sample size. The following deviations in accuracy were noted: Random Forest and mRMR models are shown to be consistently more accurate in predicting ER+, HER2-when treated with hormone therapy (both "HT" and "CT and/or HT" categories), when compared to ER-and HER2+ patients. The PAM50 basal subtype is consistently low in accuracy when testing patients treated with hormone therapy. This is most likely partially influenced by the RF and mRMR models for 'HT' to more often predict patients as sensitive, combined with the fact that ER+ and HER2-patients were more likely to response to therapy. It is important to note that the accuracy of predictions by RF and mRMR with patients treated only with chemotherapy was fairly consistent across all available subtypes (+/-10% accuracy As written, it appears to disregard the long history of predicting clinical outcome from gene expression involved in developing these classifiers from gene expression data (e.g., van't Veer et al., 2002) into clinical assays based upon expression of smaller numbers of genes.
Response: We have removed the indicated sentence, which we agree was insufficient to the comment from the previous iteration of this article: "Must be discussed in the context of existing genomics classifiers for breast cancer (e.g., OncotypeDx and/or Mammaprint)".
We in no way meant to ignore the long history of predicting clinical outcome from gene expression (as well as other genomic factors). A discussion on this topic was not included in earlier submissions as it initially had an imposed word length limit (upon first submission). We did, however, reference other articles which do discuss this topic. In Dorman (2016), which et al. described some of the methodology for initial gene selection that this study was based on, these contributions are well-referenced, including the history of the prediction of clinical outcome from genomic status: "Previous studies have derived associations between the genomic status of one or more genes and tumor response to certain therapies ( In response to Dr. Fertig's comments, we have added a short discussion with citations of previously published approaches (including MammaPrint and Oncotype DX): "Genomic information has been shown to correlate with tumor therapy response in previous studies . From these studies, analytical methods have been used to develop gene signatures for chemotherapy resistance prediction , subtypes (PAM50), and metastatic risk stratification (Oncotype DX™, MammaPrint )." Comment 5. Based on the previous review, the authors include context with other predictions of the METABRIC data in the response to the reviewers. This must also be included in the Conclusion to assess the relevance of their findings in the literature.
Response: We have added the indicated text from the previous 'response to the reviewers' (modified) to the Conclusions: "We also examined the method exhibiting the best performance in the Sage Bionetworks / DREAM Breast Cancer Prognosis Challenge , which was also phenotype-based, however it produces outcome signatures based on molecular processes, rather than the cancer drugs themselves. While interesting and informative, the results cannot be directly compared." ® 5,12-16 5 ® 1.

2.
3. outcome signatures based on molecular processes, rather than the cancer drugs themselves. While interesting and informative, the results cannot be directly compared." Please note that the majority of entries in the DREAM project were not fully curated and only exist as source code. Analyzing these files to determine what methodology was attempted by these groups is beyond the scope of our study. A description of the second place of the METABRIC phase of the DREAM challenge is provided in the link below. This link describes how the METABRIC data is trained using a bipartite graphing as input for linear models, boost models, and RankSVM. While they state that RankSVM was the least successful between the three methods, it does not appear that this particular study has been published to the literature. As a result, we cannot fully review their results, and thus cannot be compared to our methodology in the main manuscript. This study proposed prediction methods using SVM and RF classifiers with mRMR selected feature sets from cell line data and demonstrate its prediction ability for outcomes from METABRIC patient cohort. The classifiers with good prediction performance show the usefulness of combining domain knowledge with feature selection techniques. However, some details essential for reproducibility and interpretation are missing. Required information is listed in the following.
What are the values of parameters for SVM and RF classifiers and the methods for parameter selection (by default or other selection methods)?
The development and evaluation of models for patient data are not clear. Whether the models were trained using partial data from METABRIC or only leave-one-out cross-validation was applied? If cross-validation is the case, then what is the model offered at the online server because there will be more than one models created, and whether the cross-validation is involved in the feature selection process that often leads to an overestimation of the performance. For the case of training on partial data, both training and test performance are essential information for evaluating the robustness of models.
Since some of the datasets are highly imbalanced, the numbers of positives and negatives, as well as sensitivity and specificity are more important than accuracy for interpreting the results as a high accuracy with a low AUC could be the result of all positive/negative predictions on an imbalanced 3.
as sensitivity and specificity are more important than accuracy for interpreting the results as a high accuracy with a low AUC could be the result of all positive/negative predictions on an imbalanced dataset. Listing all the information along with the accuracy and AUC will help the interpretation of prediction performances.
No competing interests were disclosed.

Competing Interests:
I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above. In regards to parameter selection, the first paragraph of the methods now describes C and Sigma selection as a grid search to find the values with the lowest cross-validation misclassification rate. Similarly for RF, a grid search was used to optimize the maximum number of randomly selected genes for each tree (second paragraph of Methods section).

Comment 2:
The development and evaluation of models for patient data are not clear. Whether the models were trained using partial data from METABRIC or only leave-one-out cross-validation was applied? If cross-validation is the case, then what is the model offered at the online server because there will be more than one models created, and whether the cross-validation is involved in the feature selection process that often leads to an overestimation of the performance. For the case of training on partial data, both training and test performance are essential information for evaluating the robustness of models.
Response: We obtained new results for both RF and mRMR+SVM models when we use discovery set as training set and validation set as test set, the performance of the model was poor. After more investigation we found that there happened to be a large variation between gene expression of 26 targeted genes between discovery and validation set (please see Supplementary Dataset 2). Hence, building any classifier using discovery and validation set as training and test set in their current forms will result of poor performance, since the training and test sets are vastly different.
However, we did carry out another experiment on discovery set solely and used 70% of data for training and remaining 30% for test the performance of the model. The results have been added to the manuscript (Tables 4 and 5).
Comment 3: Since some of the datasets are highly imbalanced, the numbers of positives and negatives, as well as sensitivity and specificity are more important than accuracy for interpreting the results as a high accuracy with a low AUC could be the result of all positive/negative predictions on an imbalanced dataset. Listing all the information along with the accuracy and AUC will help the interpretation of prediction performances.
Response: As previously mentioned, we have added more performance measures including MCC and AUC. They have been added Tables 1-5 of the manuscript. 70% of the data and then tested using the remaining 30% of data as test set. We repeated this procedure 100 times and took the median as the final performance result. The results are presented in Tables 4 and 5 of the manuscript.
Comment 3: Classifiers must be put in the context of other existing genomics classifiers used in breast cancer and/or previously published in Mammaprint data.
Response: We have added two sentences in the second paragraph of the "Results and Discussion" section which describes the comparison of our gene signature to those from MammaPrint and Oncotype Dx. Pair-wise comparison of these three signatures show that they are nearly independent of one another.

Comment 4: Abbreviations SVM and RF must be spelled out as Support Vector Machine and Random Forest on first use in Methods.
Response: We thank the reviewer for this suggestion. It has been addressed in the Methods section of the manuscript.
Comment 5: Writing in SVM learning subsection of Methods requires clarification to distinguish which of these methods were developed in the previous Molecular Oncology publication and which were developed as part of this publications.
Response: This is now clarified within the first paragraph of the Methods section in the manuscript. The SVM classifier was adopted from previous Molecular Oncology publication, while the feature selection method has been developed as part of this publication.
Comment 6: Details about the SVM learning algorithm are included in the caption to Figure 1, but must also be included and completely described in text for the corresponding section of the methods.
Response: Thanks for the reviewer's suggestion. This description of the SVM learning algorithm has been moved from the Figure 1 legend and integrated into the first paragraph of the methods section.

Comment 7: No equations are provided to describe the role of the parameters C and sigma. It is also unclear whether this greedy search is implemented by the Matlab function fitcsvm or uses custom code developed by the authors.
Response: A brief description of the role of each parameter has been added to the first paragraph of the methods section of the manuscript. Readers are also now directed to a reference (Ben-Hur and Weston, 2010) if more detail is desired.
The greedy search, also called sequential backward feature selection, was implemented as a script by our lab in MATLAB. It is not a MATLAB function. This is clarified by changing a few words in the first paragraph of the methods section: "A backwards feature selection (greedy) algorithm was designed and implemented in MATLAB in which…" Moreover, as described above, the SVM classifier was adopted from previous Molecular Oncology