Gene signature‐based prediction of triple‐negative breast cancer patient response to Neoadjuvant chemotherapy

Abstract Neoadjuvant chemotherapy is the current standard of care for large, advanced, and/or inoperable tumors, including triple‐negative breast cancer. Although the clinical benefits of neoadjuvant chemotherapy have been illustrated through numerous clinical trials, more than half of the patients do not experience therapeutic benefit and needlessly suffer from side effects. Currently, no clinically applicable biomarkers are available for predicting neoadjuvant chemotherapy response in triple‐negative breast cancer; the discovery of such a predictive biomarker or marker profile is an unmet need. In this study, we introduce a generic computational framework to calculate a response‐probability score (RPS), based on patient transcriptomic profiles, to predict their response to neoadjuvant chemotherapy. We first validated this framework in ER‐positive breast cancer patients and showed that it predicted neoadjuvant chemotherapy response with equal performance to several clinically used gene signatures, including Oncotype DX and MammaPrint. Then, we applied this framework to triple‐negative breast cancer data and, for each patient, we calculated a response probability score (TNBC‐RPS). Our results indicate that the TNBC‐RPS achieved the highest accuracy for predicting neoadjuvant chemotherapy response compared to previously proposed 143 gene signatures. When combined with additional clinical factors, the TNBC‐RPS achieved a high prediction accuracy for triple‐negative breast cancer patients, which was comparable to the prediction accuracy of Oncotype DX and MammaPrint in ER‐positive patients. In conclusion, the TNBC‐RPS accurately predicts neoadjuvant chemotherapy response in triple‐negative breast cancer patients and has the potential to be clinically used to aid physicians in stratifying patients for more effective neoadjuvant chemotherapy.


| INTRODUCTION
Neo-adjuvant chemotherapy (NCT) is the standard-of-care for breast cancer and can improve treatment options for patients with large, inoperable, or advanced tumors. 1 Multiple clinical trials have illustrated the clinical benefits of NCT; in large or inoperable tumors, it has been shown to significantly reduce tumor size and enable conservative breast surgery for certain patients. The survival benefit of NCT is identical to adjuvant chemotherapy (the traditional option) in early-stage breast cancer patients. [2][3][4][5] However, in contrast with adjuvant chemotherapy, which is given in the absence of any measurable parameter for response evaluation, the response to NCT is often evaluated by MRI or PET/CT and can be classified as a pathologic complete response (pCR). Of note, patients with pCR have prolonged survival compared to patients with residual disease. 2,6,7 Despite its clinical advantages, however, only 30% of patients responded to NCT; the majority presents with residual disease (RD) and suffers from side effects that can hamper further surgical options.
To address this issue, a number of genetic signatures have been proposed to predict patient response to NCT and inform treatment decision. [8][9][10][11][12][13] In ER-positive breast cancer patients, several cell cycle pathway associated gene signatures have been commercialized due to their high accuracy in predicting response to NCT. One of the most widely utilized signatures is Oncotype DX, which predicts NCT response in ERpositive breast cancer patients based on the expression of 21 genes. [14][15][16][17][18] Patients with high Oncotype DX risk scores are more likely to respond to the NCT. 17 Other gene assays, such as EndoPredict and PROSIGNA, were also introduced to predict the NCT response in ER-positive patients. 19,20 While the success of ER-positive commercialized gene signatures has been promising, very few predictive gene signatures in ERnegative patients have been reported.
TNBC is a subset of ER-negative breast cancer, which accounts for 10-20% of all breast cancers. TNBC tumors fail to express estrogen receptors (ER), progesterone receptors (PR), and the epidermal growth factor receptor-2 (Her2). 21,22 Compared to other subtypes, TNBC is the most aggressive and is characterized by larger tumor size, higher grade, increased number of lymph node metastases at diagnosis, and the worst survival outcomes. Unfortunately, current treatment options for TNBC are very limited. 21,23,24 Indeed, no targeted therapies for TNBC are available, with the exception of PARP inhibitors in germline BRCA1/2-mutated tumors. 25 This makes chemotherapy the only treatment option for most TNBC patients. Due to the presence of a high intertumoral heterogeneity, the same NCT regimen may yield diverse responses in different patients. 26 This presents a need for the identification of predictive biomarkers that can be applied to help tailor care. Previously, Stover et al reported that both proliferation and immune-related gene signatures are associated with response in TNBC patients. 27, 28 Farmer et al reported that a stroma-related gene signature was predictive of NCT response in TNBC patients. 29 However, compared to the Oncotype DX prediction accuracy in ER-positive patients, none of these signatures achieve a high prediction accuracy, which limits their potential for clinical utilization. Currently, there is no clinically predictive gene signature for TNBC patients. 27,30 Therefore, developing clinically applicable biomarkers for TNBC to predict NCT response is critical and would spare nonresponder patients from experiencing severe side effects.
In this study, we propose a computational framework to define a whole-transcriptome signature to quantify the probability of a patient to respond to NCT. To this end, we utilized pretreated patient gene expression data by comparing the NCT responders vs nonresponders to define the gene signature associated with NCT response. Our rationale is that treatment response in cancer involves complicated cellular and molecular interactions in the tumor environment in which, for example, cell metabolism and cell-cell interactions are important. While most published signatures focus on a single tumor-associated pathway, we included all genes to capture the complicated cellular and molecular interactions, which more accurately predicts NCT response. Here, we present NCT response associated gene signatures to calculate response probability scores (RPS) in breast cancer patients.
Our results indicate the utility of our computational framework for identifying novel predictive biomarker(s) and have identified a powerful biomarker for NCT response prediction in TNBC.

| Breast cancer gene expression datasets
The Gene Expression Omnibus (GEO) and MD Anderson Cancer Center public databases were queried for available gene expression datasets using the following search terms: (breast cancer) AND (preoperative chemotherapy OR neoadjuvant chemotherapy). Only microarray datasets generated using Affymetrix (U133 and U133Plus2.0 arrays) and having more than 80 samples were included to limit the cross-platform variability. Patient samples were excluded if the biopsies were obtained after NCT, if the patient sample did not have ER-status or Her2-status information, if pathologic response was not available, or if comparable treatment agents were not found ( Figure S1). Duplicate records were removed by careful review of GEO annotations. Based on these criteria, we identified seven datasets. The GEO accession numbers and the dataset downloaded from the MD Anderson Cancer Center public database were: GSE25055, GSE20194, GSE25065, GSE20271, GSE32646, GSE22093, and Hess et al, [31][32][33][34][35][36] with sample sizes of 306, 278, 182, 178, 115, 97, and 129 respectively (Table S1 and Figure S1). In total, 115 patient expression profiles were measured by Affymetrix U133Plus2.0 arrays and 1170 patient expression profiles were measured by Affymetrix U133 array. The gene expression data were downloaded as matrices containing the expression level of all probes and then converted into gene-level expression. For genes with multiple probesets, the probeset with the largest average expression value across samples was selected to represent the expression of that gene.
GSE25055 was used as training dataset for constructing the RPS and TNBC-RPS signatures, while the other datasets were used as validation datasets. We constructed a validation metadata dataset by applying quantile normalization to re-scale the RMA normalized gene expression and then applying the ComBat function ("sva" R package) 37 to remove batch effects ( Figure S2).

| Define RPS and TNBC-RPS gene signatures
To capture the transcriptome difference between pCR and RD patients, the RPS gene signature was defined by identifying differentially expressed genes between pCR and RD patients (Table S2). A logistic regression model was constructed for each gene using patient class as the response variable (Y = 1 for pCR patients, and Y = 0 for RD patients).
Log2-transformed gene expression data were included as a predictive variable in the model (X 1 ). The coefficients (β 1 -values) and statistical significance (P-values) for each gene were estimated by applying these models to the training data (GSE25055). Then, given these values (β, P) for all genes, the RPS gene signature was constructed by using a pair of weight profiles, w + and w -, that assigned all genes which had two weights in the following way: for gene i, , where I represents the indicator function. Weights were trimmed at 10 to avoid extreme values and transformed into values within [0,1] by subtracting the minimum value and then dividing by the range. If a gene i was more significantly up-regulated in pCR vs RD samples, it received a high weight in the w + i profile and a weight of zero in the w − i profile. Conversely, a more significantly down-regulated gene in pCR vs RD samples was assigned a high weight in the w − i profile and weight of zero in the w + i profile. The TNBC-RPS gene signature was derived based on the same framework, but the logistic regression model was performed for each gene in TNBC patients only (Table S3).

| Calculation of RPS and TNBC-RPS in pretreated breast cancer samples
Given the expression profiles for a number of breast cancer patients, sample-specific RPSs were calculated for all samples based on the RPS gene signature. Specifically, a modified version of a statistical method called BASE [38][39][40][41][42] was applied as follows: first, gene expression profiles were median normalized to relative gene expression for each gene across samples. Second, for each sample, its gene expression profile was sorted in a descending order based on the relative expression to obtain an expression profile (e 1 , e 2 , …, e g ), where g was the total number of genes. The skewed distribution of up-regulated (with large values in w + ) and down-regulated (with large values in w − ) genes in pCR and RD samples were examined by comparing two cumulative functions, a foreground f(i) and a background b(i): If genes with large weights in w (w + i for up-regulated genes and w − i for down-regulated genes in breast cancer samples) also had high gene expression values in a breast cancer sample expression profile e, f (i) would accumulate more rapidly than b (i) as i increases. Third, for genes in w + i , RPS + was defined as the maximum deviation between the f (i) and b (i) and then normalized against null distribution that was generated by 1,000 iterations of a randomized tumor expression profile. The same process was applied for genes in w − i to generate the RPS -. The final RPS was determined by taking the difference between RPS + and RPS -(RPS + -RPS -). Using this approach, patients receiving high RPSs had profiles similar to gene expression profiles of patients with known pCR, while patients receiving low RPSs had profiles similar to gene expression profiles of patients with known RD.
For the TNBC-RPS calculation, the TNBC-RPS signature was applied in the TNBC patients. Following the method above, the foreground f(i) and background b(i) functions were used to calculate TNBC-RPS for each TNBC patient. Specifically, for global prediction power comparison, we calculated RPS and TNBC-RPS based on the expression of metadata. For individual cohort prediction power comparison, we calculated RPS and TNBC-RPS based on the original normalized expression data.

| Previously defined predictive gene signature calculation
Gene signatures were collected from published studies describing a variety of biological processes implicated in chemosensitivity or resistance. Three categories of gene signatures were collected in our study for comparison: Category 1 (ER-positive patient): Commercialized gene signatures were used for prediction and comparison. Oncotype DX risk scores, 8 MammaPrint signature scores, 9 EndoPredict scores, 43 Gene76 scores, 44 Genomic Grade Index (GGI), 45 and risk of recurrence scores (RORs) 46 were calculated using the "oncotypeDX", "gene70", "endoPredict", "gene76", "ggi", and "rorS" functions, respectively, from the genefu R package. 47  where w i referred to the weight of the genes in the module and e i referred to the expression of these genes. Category 2 (ER-negative and TNBC patient): The search terms: (predict OR biomarker) AND Breast cancer AND (ER negative OR Triple negative) AND (neoadjuvant OR preoperative chemotherapy) were used to find relevant publications. After excluding publications with no gene expression-based signature, or which were not validated in at least two independent datasets, 19 gene signatures remained. 28,[49][50][51] The methods of calculating those 19 gene signatures were as follows: Signature 1: Witkiewicz et al reported that cell-cycle-related genes are important for NCT and used nine genes to quantify the related pathway activity. 49 The average expression of these nine genes were calculated as the metric for prediction.
Signature 2: Turner et al presented a Consensus Signature 50 that captured the combined effect of immune function, tumor proliferation, and the tumor proliferation regulators. In detail, this signature was composed of the sum of the STAT1 module score (immune function), TOP2A (tumor proliferation), and LAPTM4B (tumor proliferation regulator) gene expression. The Module score was calculated used the equation above. We then scaled the Module score to have an inter-quartile range of 1 and a median of 0. The expression level of TOP2A and LAPTM4B was rescaled by the same method. The final score was calculated as the sum of these three scaled scores.
Signature 3: Desmedt et al 51 combined the modules associated with different tumor microenvironment components for prediction. Module scores of the Immune response, Stromal signature, and TOP2A signature (cell proliferation) were calculated through the equation described above. Specifically, the application of the signature was determined by the Her2 status. In ER-negative/Her2-negative patients, the final score was calculated as the sum of the Immune response, Module score, and Stromal Module score. In ER-negative/Her2positive patients, the final score was calculated as the sum of the Immune response, Module score, Stromal Module score, and the TOP2A signature Module score.
Signatures 4-13: Ignatiadis et al 28 reported 10 of 17 signatures that have been examined to be predictive in ER-negative patients. Similarly, the Module score of each signature was calculated through the equation above.
Signature 18: Juul et al 52 identified that the mitotic and ceramide modules were associated with the pCR and defined the paclitaxel response metagene score as the difference between mitotic Module score and ceramide Module score.
Signature 19: Farmer et al 29 used the stromal-cell-associated signature for prediction, which was calculated as the average gene expression of 48 genes.
Category 3 (Non-ER-status dependent): Stover et al 48 reported and summarized 125 signatures from previous studies for NCT prediction. For each gene signature, its Module score was calculated as the metric for prediction. In summary, a total number of 143 signatures were calculated, as described in the accompanying publications, and were validated in corresponding datasets from the original studies (Table S4-S5). Specifically, for global prediction power comparison, we calculated 143 signature scores based on the expression of metadata. For individual cohort prediction power comparison, we calculated 143 signature scores based on the original normalized expression data.

| NCT response prediction
Patients were predicted to have pCR or RD based on scores derived from the RPS, the TNBC-RPS, and the other 143 signatures collected from previous publications. For each signature, we ranked patients based on signature scores from low to high. For each patient, a threshold was set, beginning with the lowest score, where patients with a score higher than the threshold were predicted to be pCR and patients below the threshold were predicted to be RD. The sensitivity and specificity were then calculated for each threshold by comparing the predicted pCR to the actual pCR. Prediction accuracy of each signature was represented by calculating the area under the receiver operating characteristics curve (AUC).
To evaluate the performance of each signature in combination with established clinical predictors, a Random Forest model was trained to predict pCR and RD status using the RPS, the TNBC-RPS, and other signatures as predicting features, integrated with clinical predictors including age, tumor stage, and tumor grade. Random Forest classification was performed in R through the randomForest package, while setting the sample size of pCR and RD patients to be equal. 53 The performance of the model was evaluated by a 10-fold cross validation, where samples were randomly divided into 10 subgroups, with nine subgroups being used to train the Random Forest model and one subgroup used for NCT response prediction. To make each sample part of the validation set at least once, this process was repeated 10 times. Model prediction accuracy was evaluated by calculating AUC. This overall cross-validation procedure was repeated 100 times to obtain an overall average AUC.

| Pathway enrichment analysis and tumor microenvironment component decomposition
The MsigDB C2 dataset 54 was downloaded for pathway enrichment analysis. KEGG gene sets, BioCarta genes sets, and Reactome gene sets were chosen for analysis. Gene sets with less than 20 genes were excluded, which lead to the inclusion of 798 pathways. For each pathway gene set, the enrichment score was calculated based on the rank of pathway genes in the RPS and TNBC-RPS gene signatures. Specifically, the enrichment score was calculated through a walking sum method: where g i referred to the accumulative hits of genes in the gene set, d i referred to the gene rank difference between two continuous hits in the RPS or TNBC-RPS signatures, n referred to the total number of genes in the gene set, and N referred to the total number of genes in the RPS or TNBC-RPS gene signatures.
The tumor microenvironment was decomposed into three general components: infiltrating immune cells, stromal cells, and tumor cells. In detail, the abundance of infiltrating immune cells and stromal cells in the tumor microenvironment were estimated using the ESTIMATE package in R. 55 The proliferation rate of tumor cells was estimated using the normalized expression level of MKI67. 56

| Overview of the study
We developed a computational framework that could be utilized to identify predictive gene signatures associated with neoadjuvant chemotherapy (NCT) response in triple-negative breast cancer (TNBC) and then conducted a series of analyses as summarized in Figure 1. We compared pretreatment gene expression profiles between pathologic complete response (pCR) and residual disease (RD) patients from a prospective clinical study (GSE25055) to identify a weighted whole-gene signature associated with NCT response, where genes are weighted based on their capacity to discriminate pCR vs RD patients. A response probability score (RPS) was calculated for each patient in the metadata (see methods) through a rank-based algorithm called BASE. 38 Patients having high similarities between their gene expression profile and NCT response associated signature would have high RPS scores, leading to high probability of being pCR. After illustrating the efficacy of the framework by showing that the RPS has similar predictive power as the leading commercialized signatures, such as Oncotype DX and MammaPrint in ER-positive patients. We expanded the framework in TNBC patients, generated a novel TNBC response-associated signature, calculated TNBC response probability scores (TNBC-RPS) for TNBC patients in the metadata, and examined its predictive power in TNBC. Moreover, we annotated RPS and the TNBC-RPS by correlating the scores with immune cell infiltration, stromal cell abundance, and tumor cell proliferation rate in the tumor microenvironment.

| The RPS predicts patient response with high accuracy in ER-positive breast cancer
A number of gene signatures have been proposed for ERpositive breast cancer, including several commercialized assays such as Oncotype DX. [8][9][10][11][12][13] We first tested our framework in ER-positive breast cancer by comparing its performance with commercialized assays. The efficacy of our developed computational framework was validated by investigating the predictive power of the RPS for NCT response ( Figure S3A-B and Figure 2). As shown in Figure 2A, patients with pCR had a significantly higher RPS than patients with RD (P = 7e-35, Figure 2A). Because ER-negative patients had a higher response rate to NCT than ER-positive patients, 57 we separated the patients by ER status. In both ER-positive and -negative patients, pCR patients had a significantly higher RPS than RD patients (P = 5e-14, ER-positive patients; P = 9e-7, ER-negative patients; Figure 2A). Similar results were observed in the enrichment analysis, that pCR patients were significantly enriched in the high RPS group compared to other groups (P = 1e-28, All patients; P = 1e-12, ERpositive patients; P = 1e-4, ER-negative patients Figure 2B). Furthermore, to quantify the predictive power of the RPS, we utilized the RPS as a predictor to classify patients as pCR or RD. As shown in Figure 2C, RPS was predictive to NCT response and Higher predictive power was observed in ERpositive patients (AUC = 0.77, All patients; AUC = 0.78, ER-positive patients; AUC = 0.64, ER-negative patients, Figure 2C and Table S6).
To compare the predictive performance of RPS to the commercialized signatures, we collected 143 predictive signatures in breast cancer from previous publications (see methods). These signatures could be stratified based on their applicable range, including ER-positive-specific signatures, ER-negative-specific signatures, and nonspecific signatures (see methods). We applied all collected signatures in the metadata to examine and compare their AUC with the RPS in ER-positive and -negative patients. Compared with ERpositive-specific predictive signatures, the RPS had similar or higher AUC performance in ER-positive patients compared to most of the ER-positive-specific signatures ( Figure 2D and Table S6). Interestingly, MammaPrint and Oncotype DX had an AUC of 0.78 and 0.77 in ER-positive patients, respectively, while the RPS had an AUC of 0.78 in ER-positive patients ( Figure 2D and Table S6). For convenience, we grouped ERnegative-specific and nonspecific signatures together and named this group "ER-positive-nonspecific signatures". In these nonspecific signatures, the loss of RB1 expression, a cell proliferation signature, and the epithelial-mesenchymal transition (EMT) signature had the highest AUC of 0.76 in ER-positive patients. Notably, no robust signatures were identified in ER-negative patients ( Figure 2D and Table S6).
To show that the predictive accuracy of RPS in the metadata was not driven by a single dataset, we then examined and compared the predictive consistency of RPS with 143 other signatures in ER-positive patients from each dataset. As shown in Figure 2E, the RPS was predictive of the response in each individual dataset, with the lowest AUC = 0.71 in the GSE20271 dataset. This was similar to other commercially available predictive signatures, including Oncotype DX (lowest AUC = 0.69 in GSE20271) and MammaPrint (lowest AUC = 0.69 in GSE20271) (Table S6). However, the predictive ability of RPS in ER-negative patients was relatively lower compared to its prediction ability in ER-positive patients with the lowest AUC = 0.64 in the Hess et al dataset ( Figure 2F and Table S6). In summary, we validated the efficacy of our computational framework by showing the RPS's predictive power in breast cancer, particularly its comparative F I G U R E 1 The schematic diagram of this study. GSE25055 microarray data were used to determine the gene signature that captures the gene expression difference between pCR and RD patients. The signature was applied to the breast cancer metadata to calculate the patient-specific response probability score (RPS). RPS predicts response better in ER-positive patients than ER-negative patients. Following this, the TNBC gene signature was defined by using TNBC only in the GSE25055 microarray data. The signature was applied to the TNBC patients in the metadata to calculate the TNBC response probability score (TNBC-RPS) and was validated in the TNBC. The annotation of the gene signatures was performed by correlating the RPS or TNBC-RPS with the immune cell, stromal cell abundance, and tumor cell proliferation rate in the tumor microenvironment prediction power with commercialized signatures in ERpositive patients.

| The TNBC-RPS predicts NCT response in TNBC patients
After showing the efficacy of the framework in ER-positive breast cancer, we aimed to define a signature that could predict NCT response of ER-negative patients. Here, we focused on triple-negative breast cancer (TNBC), an aggressive and heterogeneous subtype. No clinically practical signatures are currently available for predicting patient response to NCT in these patients. We applied our framework to TNBC patients in the previous training dataset (GSE25055) and built a TNBCspecific signature to capture gene expression differences between pCR and RD of TNBC patients. Unsurprisingly, the TNBC-RPS calculated in the training dataset is predictive of NCT response (P = 5e-11, AUC = 0.87, Figure S3C-D). We then integrated the TNBC-RPS signature with TNBC patient expression profile in the validation metadata to calculate the TNBC-RPSs. The pCR patients had significantly higher TNBC-RPS than RD patients (P = 3e-12, Figure 3A). Moreover, the pCR patients were significantly enriched in the high-TNBC-RPS group compared to other groups, with a pCR rate of 61.2% compared to a baseline pCR rate of 33.2% (P = 5e-10, Figure 3B). We further quantified the predictive power of the TNBC-RPS in TNBC patients from our    Figure 3C). Also, the prediction power of TNBC-RPS could be observed in each individual dataset (Table S7 and S8). The performance of the TNBC-RPS to predict NCT response was compared to previously defined predictive signatures. As stated in the previous section, we collected 143 predictive signatures, which were composed of ER-positivespecific signatures, ER-negative-specific signatures, and nonspecific signatures. From these, 19 ER-negative-specific signatures were identified, and the other 124 signatures were grouped into ER-negative-non-specific signatures for comparison. As shown in the upper panel of Figure 3D, the TNBC-RPS outperformed 19 ER-negative-specific predictive signatures in predicting NCT response, with an AUC of 0.77. The next-highest AUC was achieved by the loss of PTEN gene signature (AUC = 0.67, Table S7), which has been reported to predict NCT response in TNBC patients. 28,58 The TNBC-RPS also outperformed all ER-negative-nonspecific predictive signatures in the validation metadata ( Figure 3D and Table S7). The highest AUC of the ER-negative-nonspecific predictive signature was achieved by an E2F1 pathway-related gene signature (AUC = 0.68, Table S7). We further examined the AUC of each signature across the individual dataset ( Figure 3E) and found that the TNBC-RPS was predictive of NCT response in TNBC patients across all datasets, with the highest AUC = 0.91 in GSE20194 and Hess et al dataset (Table S7). In three of six datasets, the TNBC-RPS had an AUC higher than 0.75 while other signatures did not present such consistent prediction power ( Figure 3E and Table S7). To compare the predictive accuracy between the TNBC-RPS and other signatures, a paired T-test was used to measure the statistical significance of AUC differences across six individual datasets. As shown in Figure 3F, the TNBC-RPS significantly outperformed 129 of 143 signatures in predicting NCT response (P < .05, Figure 3F). In cases in which the TNBC-RPS was not statistically significant compared to other signatures, a positive trend in the T-score was observed; this still indicates a better prediction ability when using the TNBC-RPS as the predictor.

| The TNBC-RPS predicts NCT response in each clinical stage and grade
Although the TNBC-RPS showed good prediction power in TNBC patients, we were concerned that tumor stage or grade might have confounded these findings; it has been reported that TNBC patients with a more advanced tumor stage or grade tend to have better response to NCT. 59 To evaluate this, we first examined the predictive ability of the TNBC-RPS in the validation metadata for each tumor stage. By calculating the TNBC-RPS for each individual stage in both pCR and RD patients (Table S9), we found that pCR patients had significantly higher RPS than RD patients (P = .007, Stage I; P = 9e-6, Stage II; P = .004, Stage III; P = 1e-6, Stage IV; Table S9). Secondly, we calculated the AUC of the TNBC-RPS in a stage-specific manner. As shown in Figure 4, the TNBC-RPS could predict stage-specific responses in TNBC patients, indicating that the predictive power was not affected by stage stratification (Figure 4A-D). Interestingly, the TNBC-RPS showed a high predictive power within stage-I patients (AUC = 0.82, TNBC-RPS; Figure 4A), indicating that the TNBC-RPS could robustly predict NCT response in early-stage breast cancer patients. This is important since the development of diagnostic techniques increases the number of patients diagnosed at early stages, thus requiring predictive markers that are effect at those stages.
Similar to tumor stage, tumor grade may also confound the prediction of the TNBC-RPS in TNBC patients. Hence, we performed the same analyses by using the TNBC-RPS as the predictor in each tumor grade and found that AUC = 0.64 and AUC = 0.75 in grade-II and grade-III patients ( Figure 4E-F) respectively. The conclusion of the TNBC-RPS being a grade-specific predictor was further validated by comparing the TNBC-RPS difference between pCR and RD patients (Table S9).

| The TNBC-PRS adds additional predictive power to current clinical predictors
We have demonstrated that the TNBC-RPS could predict pCR in each TNBC clinical stage and grade. In clinical practice, the combination of clinical stage, grade, and age is used to predict NCT response. 60 Therefore, we investigated whether adding the TNBC-RPS to current clinical predictors could further improve prediction accuracy. First, we applied Random Forest algorithm to calculate AUC for clinical predictors in TNBC patients and performed 10-fold cross-validation. As shown in Figure 5A, the prediction accuracy only achieved an AUC of 0.69 with the use of clinical predictors.
Then, by adding the TNBC-RPS to clinical predictors, we improved the AUC to 0.80 (Table S10). In order to further understand the contribution of the TNBC-RPS to the prediction, we investigated the relative importance of the TNBC-RPS and clinical factors through a 10-fold cross validation process. As expected, the relative importance of the TNBC-RPS was significantly higher than other clinical predictors, which indicated the predominant predictive power of the TNBC-RPS in the TNBC NCT response prediction (P < 2e-16, Figure 5B).
Next, we assessed whether other gene expression signatures displayed similar properties. Thus, we repeated the same analysis as mentioned above and combined clinical predictors with each of the 143 signatures to test the AUC change. As shown in Figure 5C, 90 of 143 signatures had an AUC higher than 0.7, with TNBC-RPS having the highest AUC = 0.80. The second-highest signature was reported by Witklewicz et al, with an AUC = 0.74 (Table S10). In summary, the TNBC-RPS combined with the clinical predictors outperformed the prediction accuracy compared to the previous signatures (AUC = 0.80, Figure 5C and Table S10). In addition, we also performed the same analyses by using RPS in ER-positive patients and found that RPS could further improve the prediction accuracy of the current clinical predictors to AUC = 0.79 (Table S11).

| The TNBC-RPS associates with immune cell infiltration, stromal cell abundance, and cell proliferation
To biologically annotated RPS and TNBC-RPS, which could potentially indicate the different mechanisms underlying the chemotherapeutic response between ER-positive and TNBC patients, we performed a pathway enrichment analysis in both gene signatures ( Figure 6A and Table S12). A positive enrichment score refers to the enrichment of pathways in the up-regulated genes of the signature, while a negative enrichment score refers to the enrichment of pathways in the down-regulated genes of signature. Interestingly, we found that some pathways were shared in both signatures, while some pathways were presented in a signature-specific way. Cell-cycle-related pathways were shared between the RPS and TNBC-RPS, indicating the involvement of cell-cycle pathways in the NCT response. For example, the KEGG cell cycle pathway was shared by the gene signatures of the RPS and TNBC-PRS, with an enrichment score 0.20 and 0.16 respectively (Table S12). Moreover, pathways related to immune response were also found in both the RPS and TNBC-RPS. The REACTOME antigen-presenting pathways were enriched in the TNBC-RPS gene signature (with enrichment scores of −0.12) and the KEGG T-cell-receptorrelated pathways were enriched in the RPS gene signature (with enrichment scores of −0.11). In addition to shared pathways, we also identified signature-specific pathways. For example, protein transportation-related pathways, like the REACTOME Extracellular Matrix (ECM) pathway, were only enriched in the TNBC-RPS gene signature (with an enrichment score of 0.11) (Table S12).
We next performed clustering analysis on the pathway-enrichment scores in the RPS and TNBC-RPS to better compare these two signatures ( Figure 6A). The enriched cell-cycle and immune-related pathways in both the RPS and TNBC-RPS gene signatures clustered together, while the enriched ECM pathways only formed a unique cluster in the TNBC-RPS gene signature. Moreover, by performing gene ontology (GO) enrichment analysis in the RPS and TNBC-RPS gene signatures, we validated biological processes that had been identified by the pathway analyses (Tables S13 and S14). From both the pathway and GO enrichment analyses, we hypothesized that the RPS and TNBC-RPS could capture tumor microenvironment characteristics, in which the RPS reflects both the cell-cycle and immune-related pathway activities, while the TNBC-RPS reflects the activities of pathways including the cell-cycle, immune-related, and ECM pathways.
To examine our hypothesis, we deconvoluted the tumor microenvironment into three general components-infiltrating immune cell abundance, stromal cell abundance, and tumor cell proliferation rate-to represent the activation of immune-related, ECM-related, and cell-proliferation-related pathways (see methods) respectively. The association between the RPS and those three components in ER-positive patients was examined by Spearman correlation. The RPS was positively correlated with immune cell infiltration and tumor cell proliferation rate (SCC = 0.27, Immune related A immune cell infiltration; SCC = 0.28, tumor cell proliferation rate, Figure 6B), but was not associated with stromal cell abundance (SCC = 0.03, Figure 6B). Moreover, we performed similar analyses using the TNBC-RPS and demonstrated that it was strongly negatively correlated with stromal cell abundance and was weakly correlated with immune cell infiltration and tumor cell proliferation rate (SCC = −0.57, stromal cell abundance; SCC = −0.17, immune cell infiltration; SCC = 0.11, tumor cell proliferation rate; Figure 6C). Given the fact that the TNBC-RPS mainly correlated with stromal cell abundance, we examined if stromal cell abundance could be used for prediction. We found that stromal cell abundance was predictive for the NCT in TNBC patients, with an AUC = 0.55 ( Figure S3E-F).

| DISCUSSION
Neoadjuvant chemotherapy is being used more and more frequently for treating breast cancer patients. This is due to its advantages in reducing tumor size, improving surgical options, and significantly increasing survival in responders. However, broad clinical application remains questionable because of a low response rate and the potential for significant side effects. The most extreme case is TNBC, which is the most aggressive subtype of breast cancer and has the worst prognostic outcome. Due to its heterogeneity, patients with TNBC respond differently to NCT. Numerous efforts have been put into developing the predictive signatures for TNBC, but, currently, there is no clinically applied predictive signature. Therefore, there is an urgent need for developing robust predictive biomarkers for TNBC patients. Although many studies have focused on the chemotherapy regulatory program difference between pCR and RD, 61-63 the mechanisms underlying the survival of resistant tumor cells remain poorly understood.
In this study, we developed a novel framework for identifying predictive gene signatures in breast cancer patients. We validated the efficacy of this framework by showing that the RPS predicted NCT response in breast cancer patients, particularly in ER-positive patients (Figure 2A-C and  Table S6). In addition, compared to the commercialized signatures, the RPS had a comparable prediction ability across each individual dataset ( Figure 2D-E and Table S6). We then applied the framework to TNBC patients and calculated the TNBC-RPS. The TNBC-RPS was predictive of the response in TNBC patients ( Figure 3A-C and Table S7). Compared to the previously-developed ER-negative-specific and nonspecific prediction signatures, the TNBC-RPS outperformed 143 predictive gene signatures and presented robust prediction accuracy ( Figure 3D-F and Table S7). Of importance, the TNBC-RPS leads to a higher AUC of up to 0.80 in TNBC patients ( Figure 5A-B) and exceeded the performance of the 143 predictive gene signatures when combined with clinical predictors ( Figure 5C). We, therefore, provide a new framework for identifying predictive markers of NCT response. In addition, to facilitate the clinical utility of RPS and TNBC-RPS signatures, we also provided a revised version of those two gene signatures with fewer genes (Table S15).
Previous studies calculated the scores of different gene signatures using only a single method. This strategy does not take into account the variation in the methods used to calculate the scores from the gene signatures. Instead of using this one-method-fits-all strategy, we validated the previously published signatures by applying the same algorithms that were used to calculate the scores of each signature to the same datasets and reproduced the published prediction performances. Then, we applied the gene signatures to the validation metadata for prediction. This made the prediction accuracy comparison more objective since we took the impact of different methods into consideration (Figures 2 and  3). In Table S4, we present the validation results of the previous signatures. We acquired consistent results by repeating previous published gene signatures in their validation datasets. Despite the subtle differences between the P-value reported previously and our calculated AUC (likely caused by the update or different normalization methods on the raw microarray data), we showed that our model significantly outperformed most of the reported signatures.
The drug-response mechanisms in breast cancer have been studied for many years but were still poorly understood. We investigated the association between the RPS and characteristics of the ER-positive tumor microenvironment, as well as between the TNBC-RPS and characteristics of the TNBC tumor microenvironment respectively ( Figure 6). Of note, the RPS identified changes in the tumor cell proliferation rate and immune cell infiltration in ER-positive patients, which was supported by previous studies showing that the cell-cycle-related 16,64-66 and immune-infiltration-related gene signatures [67][68][69] were associated with responsiveness. This observation could be further validated through the prediction performance of the 143 predictive gene signatures. For example, Oncotype DX, a signature composed of cell-cycle-related genes, and the Immune Signature Gene Module score were both predictive to the response in ER-positive patients (Table S6). 16,28 The TNBC-RPS primarily captured the relative abundance of the stromal cells in the tumor microenvironment. Farmer et al reported the similar finding in TNBC patients, as well. 29 Meanwhile, we also used the stromal cell abundance for prediction in TNBC patients and got an AUC = 0.55, indicating a predictive role of stromal cells in TNBC patients' NCT response ( Figure S3E-F). [69][70][71][72] Therefore, our findings provide an understanding of cancer biology in breast cancer by showing which aspect(s) of the tumor microenvironment might influence the response to the NCT.

ZHAO et Al.
Although we have demonstrated the efficacy of the RPS and the TNBC-RPS in predicting the response to NCT, the prediction power and the applicable range of the model could be further improved. In addition to the gene signatures, other IHC-staining signatures or MRI imaging-based prediction models were used to predict the response to NCT. [73][74][75] However, we lack the data to compare the performance of our signatures to these methods or to integrate them into the model for better prediction. Moreover, our signatures were applicable to the prediction of the combination of antimetabolite-, anthracycline-, alkylating agent-, and taxane-based chemotherapy-treated patients and have not been extended to investigate its predictive power with other chemotherapy agents or targeted therapy agents. With the release of more gene expression data, it may be possible to extend the applicable range of our signatures or to develop drug-specific-predictive gene signatures.
In summary, we developed a framework for identifying a predictive gene signature in breast cancer and defined two gene signatures that could be used to predict NCT response in ER-positive and TNBC patients respectively. We have demonstrated that the RPS performed at a comparable level to the current commercialized signatures, while the TNBC-RPS outperformed 143 gene signatures for TNBC patients in prediction. More importantly, integrating the RPS or TNBC-RPS with current established clinical predictors enhanced the predictive power, compared to using the clinical predictors only. In addition, the RPS and TNBC-RPS captured different aspects of the tumor microenvironment, leading to tantalizing insights as to the potential biological mechanisms driving differences in the chemotherapeutic response. This computational framework can also be readily extended to define predictive biomarkers in other cancer types.

CONFLICTS OF INTERESTS
The authors declare that they have no competing interests.

AUTHOR CONTRIBUTION
Conception and design: YZ and CC. Development of methodology: YZ and CC. Acquisition of data: YZ and CC. Analysis and interpretation of data: YZ, ES, and CC. Writing, review, and/or revision of the manuscript: YZ, ES, and CC. All authors read approved the final manuscript.