Tamoxifen therapy benefit predictive signature coupled with prognostic signature of post-operative recurrent risk for early stage ER+ breast cancer.

Two types of prognostic signatures for predicting recurrent risk of ER+ breast cancer patients have been developed: one type for patients accepting surgery only and another type for patients receiving post-operative tamoxifen therapy. However, the first type of signature cannot distinguish high-risk patients who cannot benefit from tamoxifen therapy, while the second type of signature cannot identify patients who will be at low risk of recurrence even if they accept surgery only. In this study, we proposed to develop two coupled signatures to solve these problems based on within-sample relative expression orderings (REOs) of gene pairs. Firstly, we identified a prognostic signature of post-operative recurrent risk using 544 samples of ER+ breast cancer patients accepting surgery only. Then, applying this drug-free signature to 840 samples of patients receiving post-operative tamoxifen therapy, we recognized 553 samples of patients who would have been at high risk of recurrence if they had accepted surgery only and used these samples to develop a tamoxifen therapy benefit predictive signature. The two coupled signatures were validated in independent data. The signatures developed in this study are robust against experimental batch effects and applicable at the individual levels, which can facilitate the clinical decision of tamoxifen therapy.


INTRODUCTION
Breast cancer is the most prevalent cancer among women worldwide and approximately 70% of cases express estrogen receptor [1,2]. Tamoxifen has been the major adjuvant therapy for ER+ breast cancer, but onethird of early-stage patients treated with tamoxifen after surgery for five years will experience a relapse of cancer within fifteen years [3,4]. To reduce the recurrent rate, the majority of early-stage ER+ patients also receive adjuvant chemotherapy after surgery, of which only a small proportion will ultimately benefit from the adjuvant chemotherapy, while all remain at risk of toxic side-effects [5]. Therefore, a signature for identifying patients who can benefit from tamoxifen therapy is required. In addition, although continuing tamoxifen therapy has been found to produce a reduction in recurrence and mortality for ER+ breast cancer [6], the patients treated with tamoxifen for a long time may suffer from side-effects, such as deep-vein thrombosis, endometrial cancer, pulmonary embolus, bone loss, stroke and genito-urinary system dysfunction [7][8][9]. If patients at low risk after surgery can be discriminated from patients at low risk with the help of tamoxifen therapy, clinicians could make proper decisions on tamoxifen therapy for the two distinct groups to assure its effectiveness and minimize adverse treatment effects.
Many prognostic signatures, such as the 70-gene signature reported by van't Veer et al. [10] and the 76-gene signature reported by Wang et al. [11], have been developed for predicting clinical outcome of ER+ breast cancer patients accepting surgery only [12,13]. Although these drug-free prognostic signatures could be used to guide the recommendation of adjuvant tamoxifen therapy based on the finding that only patients in the high-risk group may benefit Oncotarget 44594 www.impactjournals.com/oncotarget from tamoxifen therapy [14], they cannot further distinguish high-risk patients who cannot benefit from tamoxifen therapy. Some other researchers used samples of ER+ breast cancer patients receiving post-operative tamoxifen therapy to develop signatures for predicting clinical outcome of these patients [15][16][17][18]. Patients with low risk of recurrence recognized by such signatures are considered to be able to benefit from tamoxifen therapy and might be recommended to tamoxifen therapy. However, some of these patients will be at low risk of recurrence if they accept surgery only and actually need no tamoxifen therapy after surgery. Obviously, these problems need to be solved.
Most of previously reported signatures are based on risk scores, usually calculated as some summaries of expression measurements of the signature genes, to allocate patients into different prognostic groups [10][11][12][16][17][18]. However, such risk-score based signatures often fail in independent samples [19][20][21][22] because risk scores summarized from expression measurements of signature genes are sensitive to experimental batch effects [22][23][24]. As the applications of such risk-score based signatures require data normalization using a set of samples [10][11][12][16][17][18], the risk classification of a sample depends on the risk composition of the samples analyzed together with this sample [25]. In contrast, the signatures based on the within-sample relative expression orderings (REOs) of gene pairs are insensitive to experimental batch effects and invariable to monotonic data normalization [22][23][24][25]. Based on this unique advantage, the REO-type prognostic signatures can perform robustly in interlaboratory datasets and allow application at the individual levels [15,26]. Another important advantage of REOs is that we can pool samples from different small datasets together for further analysis, which is of special interest given that the discovery and validation of prognostic signatures often need a large number of samples [26]. However, one major problem of finding REO-type prognostic signatures is that the number of gene pairs constituted by all genes in a dataset is extremely large, leading to a super-high dimensional problem and consequently a over-fitting problem [27]. To improve the robustness of analytical results, a commonly approach is to start with pathway analyses to develop a signature based on the phenomenon that signatures identified from different samples are often closely related in functions [28,29]. Our previous research has found that within-sample REOs are overall stable in particular types of normal human tissue but widely disturbed in the corresponding cancers, which could provide the basis for pathway analysis based on REOs [30].

Drug-free prognostic signature of post-operative recurrent risk
Using the gene expression profiles of 167 normal breast tissue samples measured by the GPL96 platform (Affymetrix HG-U133A) ( Table 1), we identified 22,717,681 stable gene pairs, each of which had a stable REO in more than 99% of normal samples. Similarly, we identified 45,603,713 stable gene pairs in 407 normal breast tissue samples (Table 1) measured by the GPL570 platform (Affymetrix HG-U133 plus 2.0). The two lists of stable gene pairs had 17,507,393 overlaps, of which more than 98% had identical REOs, which was highly unlikely to occur by chance (p < 1.0E-16, binomial distribution test, see Methods). The highly stable REOs reflect the coordinated structure of gene expressions in the normal breast tissue, based on which we could characterize every cancer sample by identifying gene pairs with reversal REOs in this sample [30]. In the following text, we used the gene pairs with stable normal REOs consistently detected by both the GPL96 and GPL570 platforms to characterize cancer samples.
The 219 samples of lymph-node-negative patients accepting surgery only, collected from the GSE7390 and GSE6532_ut datasets (Table 1), were used as the discovery cohort to develop a drug-free prognostic signature of post-operative recurrent risk. Firstly, based on the 1320 canonical pathways documented in the C2 collection of the MSigDB, we identified pathways whose disrupted REOs were significantly associated with recurrence-free survival (RFS). Here, RFS was used in a broad sense to represent the prognostic end points of both local recurrence and distant recurrence [47]. For each pathway, among the intrapathway gene pairs with stable REOs in normal tissue, the frequency of gene pairs with reversal REOs in each cancer sample was calculated, termed as the disruption index of this pathway in this sample. Then, using the univariate Cox proportional-hazard model, with FDR < 5%, we identified 37 pathways whose disruption indexes were significantly correlated with RFS (Supplementary Table 1). To search for significantly correlated RFS-relevant pathways, we evaluated the correlations of the disruption indexes among the RFS-relevant pathways using Spearman rank correlation test with FDR < 5%. After linking every two significantly correlated pathways whose Spearman rank correlation coefficient was larger than 0.6, we found 23 pathways that could be connected together as a large network (Supplementary Figure 1). Many of these 23 pathways are well-known metastasisassociated pathways, including P53 and RAS signaling pathways, cell-cycle-related pathways and immunityrelated pathways, as described in Supplementary Table 1. Finally, we searched for prognostic signature of gene pairs within these 23 RFS-relevant pathways, which could be regarded as the core drug-free RFS-relevant pathways. By this way, the number of gene pairs to be searched was greatly reduced, which was expected to be able to improve the robustness of signature selection.
Within the 23 pathways, there were 19,844 gene pairs with stable REOs in the normal breast tissue. From these gene pairs, using the univariate Cox proportional-hazard www.impactjournals.com/oncotarget model, with FDR<10%, we indentified 138 gene pairs whose reversal REOs were significantly correlated with poor RFS (see Methods). From these 138 gene pairs, a forwardstepwise selection algorithm was performed to obtain a subset of gene pairs whose C-index reached maximum (see Methods) based on the following classification rule: patients with no reversal gene pairs in the subset were assigned to the low-risk group and all the other patients were assigned to the high-risk group. Finally, we extracted nine gene pairs (Table 3), termed as the drug-free prognostic signature of post-operative recurrent risk, which classified the discovery cohort into a low-risk group with 110 patients and a highrisk group with 109 patients. As shown in Figure 1A, the patients in the low-risk group had significantly better RFS than the patients in the high-risk group (HR = 3.99, 95%CI:2.47-6.45, p = 1.02E-09, C-index = 0.69).
In the first independent validation cohort of the GSE2034 dataset, the drug-free prognostic signature identified 112 patients at low risk and 97 patients at high risk, respectively, while the RFS of the former was significantly better than that of the latter (HR = 1.95, 95%CI:1.25-3.04, p = 2.71E-03, C-index = 0.59, Figure 1B). The drug-free prognostic signature was also validated in another independent GSE4922_ut dataset: the low-risk group of 82 patients had a significantly better RFS than the high-risk group of 34 patients (HR = 2.61, 95%CI:1.31-5.19, p = 4.49E-03, C-index = 0.60, Figure 1C). The first validation cohort lacks clinical data, while multivariate Cox analyses for the discovery cohort and the second validation cohort both showed the drug-free prognostic signature was a strong independent factor for predicting the postoperative recurrent risk after adjusting age, tumor size and histology grade (Table 4).
Taken together, the above results demonstrated that the drug-free prognostic signature could robustly predict recurrent risk of ER+ breast cancer patients after surgery.

Tamoxifen therapy benefit predictive signature
For samples of ER+ breast cancer patients receiving post-operative tamoxifen therapy, we firstly used the drugfree prognostic signature to recognize patients who would have been at low risk of recurrence if they had accepted surgery only, and then used the remained high-risk samples to develop a therapy benefit predictive signature for identifying patients who could benefit from tamoxifen therapy ( Figure 2).
Notably, the datasets of patients receiving postoperative tamoxifen therapy also included samples of lymph-node-positive patients ( Table 2). Under the assumption that both lymph-node-positive and lymph-nodenegative patients with high risk of recurrence after surgery would be the same likely to have micro-distant-metastases, we pooled high-risk patients predicted from both lymphnode-positive and lymph-node-negative patients together as the discovery cohort. For each of the four datasets including both lymph-node negative and positive samples, we found no differentially expressed genes (DEGs) between the highrisk patients of the lymph-node positive and negative group using Student's t-test, with FDR < 5%. Similarly, no DEGs were found between the low-risk patients of the lymphnode positive and negative group. On the other hand, we found that DEGs between the high-and low-risk groups for lymph-node negative patients was consistent with the corresponding DEGs for lymph-node positive patients. From the GSE17705 dataset, we detected 7075 and 6221 DEGs between the low-and high-risk groups for the lymphnode negative and positive patients, respectively. The two lists of DEGs shared 5312 genes and they all showed the same deregulation directions (up-or down-regulation) in the high-risk patients compared with the low-risk patients, which was highly unlikely to occur by chance ( p < 1.0E-16, binomial distribution test). Similarly, for the GSE6532_tt1, GSE6532_tt2 and GSE4922_tt datasets, the DEGs between the distinct prognostic groups for lymph-node negative Oncotarget 44596 www.impactjournals.com/oncotarget patients were also highly consistent with the corresponding DEGs for lymph-node positive patients (all p < 1.0E-16, binomial distribution test). These results provided evidence that the drug-free prognostic signature was independent of the lymph-node status.
Applying the drug-free prognostic signature to a total 521 samples of ER+ breast cancer patients receiving post-operative tamoxifen therapy, collected in the GSE17705, GSE12093 and GSE6532_tt1 datasets, we recognized a total 320 high-risk patients (184, 68 and 68 in the three datasets, respectively). These 320 patients who would have been at high-risk of recurrence if they had accepted surgery only were used as discovery cohort to develop a tamoxifen therapy predictive signature. Then, we developed the tamoxifen therapy benefit predictive signature in the same way of developing the drug-free prognostic signature. Briefly, we firstly identified 89 RFSrelevant pathways using the univariate Cox proportionalhazard model with FDR < 5% (Supplementary Table 2), and from which we further extracted 46 strongly correlated pathways that could be connected together as a large network by linking every two significantly correlated pathways whose Spearman rank correlation coefficient was larger than 0.6 (Supplementary Figure 2). We defined the 46 pathways as the core tamoxifen-associated RFSrelevant pathways and some references suggesting their relevance to tamoxifen resistance were listed in Supplementary Table 2.  [44] GPL96 85 GSE2034 [11] GPL96 209 GSE4922_ut c [45] GPL96 116 Samples of patients receiving post-operative tamoxifen therapy GPL96 136 GSE6532_tt1 b [44] GPL570 87 GSE6532_tt2 b [44] GPL96 176 GSE4922_tt c [45] GPL96 66 GSE9195 [46] GPL570 77 Note: a To determine stable gene pairs in normal tissue, from each dataset only normal samples were collected and the information of disease samples was not presented. b GSE6532 series contains three type samples: GSE6532_ut, samples of the lymph-node-negative patients accepting surgery alone; GSE6532_tt1, samples of the patients receiving post-operative tamoxifen therapy measured by GPL570 platform; GSE6532_tt2, samples of the patients receiving post-operative tamoxifen therapy measured by GPL96 platform. c GSE4922 series contains two type samples: GSE4922_ut, samples of the lymph-nodenegative patients accepting surgery alone; GSE4922_tt, samples of the patients receiving post-operative tamoxifen therapy.
The datasets in bold were discovery cohort.
Oncotarget 44597 www.impactjournals.com/oncotarget From 10,096 gene pairs with stable REOs within these 46 pathways in the normal tissue, we identified 67 gene pairs whose reversal REOs were significantly correlated with poor RFS using the univariate Cox proportional-hazard model with FDR < 10%. From these 67 gene pairs, we performed a forward-stepwise selection algorithm to extract a subset of gene pairs with the highest C-index based on the following classification rule: patients were assigned to the tamoxifen benefit group if no gene pair in the subset was reversed and all the other were assigned to the tamoxifen non-benefit group. Finally, a tamoxifen therapy benefit predictive signature consisting of ten gene pairs (Table 5) was identified, which allocated the 320 drug-free high-risk patients into a tamoxifen benefit group of 168 patients and a tamoxifen non-benefit group of 152 patients, respectively. The RFS of the former was significantly better than that of the latter (HR = 5.27, 95%CI:3.13-8.87, p = 3.03E-12, C-index = 0.70 Figure 3A).
In the first independent validation dataset GSE6532_ tt2, for the 127 high-risk patients recognized by the drug-free prognostic signature, 55 and 72 patients were classified into tamoxifen benefit and non-benefit groups, respectively, and the former had a significantly different RFS from the latter (HR = 2.99, 95%CI:1.54-5.82, p = 7.26E-04, C-index = 0.64 Figure 3B). From the independent GSE4922_tt and GSE9195 datasets, 34 and 72 drug-free high-risk patients were recognized by the drug-free prognostic signature, respectively, and we pooled them together as the second validation cohort. The therapy benefit predictive signature could stratify this validation cohort into a tamoxifen benefit group of 85 patients and a tamoxifen non-benefit group of 21 patients with significantly different RFS (HR = 3.38, 95%CI:1.65-6.92, p = 4.15E-04, C-index = 0.63 Figure 3C). In addition, for each of the discovery and validation cohorts, the RFS of the tamoxifen benefit group was not significantly different from that of the drug-free low-risk group recognized by the drug-free prognostic signature, while the latter group also had significantly better RFS than the tamoxifen non-benefit group (Figure 4). The similar results were observed when applying the two coupled signatures to lymph-node-negative and lymphnode-positive patients separately (Supplementary Table 3). A multivariate analysis in the discovery cohort was not performed due to a number of missing values, while multivariate Cox analyses for the two validation cohorts both showed that the therapy benefit predictive signature remained significantly associated with RFS after adjusting for clinical factors of age, node status, tumor size and histology grade ( Table 6).
The GSE6532 series included 85 samples of lymph-node-negative patients accepting surgery only (GSE6532_ut) and 114 lymph-node-negative patients treated with tamoxifen (GSE6532_tt1 and GSE6532_ tt2). Thus, we could compared RFS between the tamoxifen-treated and the tamoxifen-untreated patients in each of the three groups classified by the two coupled signatures. As expected, in the drug-free low-risk group, RFS of the 28 tamoxifen-treated patients were not significantly different from that of the 43 tamoxifenuntreated patients (HR = 1.14, 95%CI:0.37-3.55, p = 0.8179, Figure 5A). Also, in the tamoxifen nonbenefit group, the 49 tamoxifen-treated patients had no significant better RFS than the 21 tamoxifen-untreated patients (HR = 0.86, 95%CI:0.41 amoxifen-untreated   Figure 5C). These results suggested that both the drug-free low-risk patients and the tamoxifen non-benefit patients could not benefit from tamoxifen therapy. In contrast, in the tamoxifen benefit group, the 37 tamoxifen-treated patients had a significant better RFS than the 21 tamoxifen-untreated patients (HR = 0.41, 95%CI:0.17-0.99, p = 0.0415, Figure 5B). Similar comparison results were found in a merged dataset that included 233 samples of lymph-nodenegative patients receiving post-operative tamoxifen therapy (GSE17705, GSE4922_tt and GSE9195) and 459 samples of lymph-node-negative patients accepting surgery only (GSE2034, GSE7390 and GSE4922_ut) ( Figure 5D, 5E, 5F). This comparison analysis could not be performed for lymph-node-positive patients because there were no samples of lymph-node-positive patients without accepting tamoxifen therapy.
Taken together, the above results suggested that the two coupled signatures could be used to facilitate the clinical decision of tamoxifen therapy.

DISCUSSION
In this study, we identified a therapy benefit predictive signature coupled with a drug-free prognostic signature for early stage ER+ breast cancer patients. The two signatures can be used sequentially to stratify early stage ER+ breast cancer patients into three groups. The  first group includes patients who will be at low-risk of recurrence if they accept surgery only, and we could recommend them to accept no or a short duration of tamoxifen treatment. The second group includes patients who will be at high risk of post-operative recurrence but can benefit from tamoxifen therapy. For these patients, the decreased risk after tamoxifen therapy could be attributed to the tamoxifen efficacy, and thus tamoxifen therapy could be recommended to them. For the third group of patients who will keep at high risk after tamoxifen therapy, we can infer that the routine clinical tamoxifen therapy cannot improve their clinical outcomes. Different from previously reported prognostic signature, the two coupled signatures can find most of patients who could benefit from tamoxifen therapy and the patients at low risk with surgery only, and thus insulating them from cytotoxic chemotherapy or even tamoxifen therapy.
Notably, for the third group of patients, we should not simply infer that they are resistant to (or cannot respond to) tamoxifen. Some of these patients, who Oncotarget 44601 www.impactjournals.com/oncotarget could have poor prognoses on account of their resistance to drug-induced tumor cell apoptosis [48], could be considered to be truly resistant to tamoxifen, so prescription of other treatment modalities such as chemotherapies or target therapies could be recommended [49,50]. However, a large portion of these patients could respond to tamoxifen but the therapy efficacy may be insufficient in competition with tumor growth ability [51,52]. If this is the case, a larger dosage and longer duration of tamoxifen therapy could be recommended [53]. Thus, the therapy benefit predictive signature can be regarded as an apparently resistant signature which can be used to predict whether the prognosis of a patient can be improved by the routine clinical tamoxifen therapy.   Gene A has a higher expression level than Gene B in normal breast tissues. To identify a drug resistant signature for discriminating patients who can respond to tamoxifen, we need gene expression data of responders an non-responders of patients accepting tamoxifen therapy, which, however, are currently unavailable for post-operative patients. Nevertheless, samples of metastatic patients accepting tamoxifen therapy, whose response to the treatment can be clearly defined [54], could be subjected to gene expression profiling to develop the drug resistant signature .
In clinical practice, almost all lymph-node positive patients undergo lymphadenectomy [55] and after that they should have low risk of recurrence if they have no micro-distant-metastases. We assumed that high-risk patients predicted from either the lymph-node negative or positive group by the drug-free prognostic signature would be the same likely to have micro-distant-metastases. Thus, the signature should be independent of the lymphnode status, as evidenced by the observation that the transcriptome difference between the distinct prognostic groups for lymph-node negative samples was consistent with the corresponding difference for lymph-node positive samples and no transcriptome difference could be observed between the same prognostic groups predicted from the lymph-node positive and negative patients. All of these suggested that high-risk patients of the lymph-node positive and negative group possess similar molecular characteristics.
For clinical application, we can develop a custom array or RT-PCR kit to measure expression intensities of the 32 genes included in the two coupled signatures to determine the REOs of the signature gene pairs. Compared with the microarray technique, the RT-PCR technique is more reliable and reproducible for quantitation of transcriptional abundance of genes. Notably, the problem of experimental batch effect and data normalization also exists when RT-PCR is used to measure gene expression intensities [56]. However, it can be expected that REOs deduced from gene intensities measured by RT-PCR tend to be robust against experimental batch effects.
Due to the high-dimension problem inherent in microarray data, especially when we focus on analyzing a huge number of gene pairs, the identification of disease signatures is liable to false discoveries [27]. Through mapping gene pairs into pathways, we started with pathways to improve the robustness of the identification of signatures. As demonstrated in this study, the identified signatures can perform robustly in independent datasets. However, due to the limited gene annotation to biological pathways [57,58], some important pathways associated with survival might be missed. A method worth exploring is to augment annotated genes of pathways using genes that are closely linked with intra-pathway genes in proteinprotein interaction network [59,60].
In this study, in order to ensure the robustness of signature performance in samples detected by different Affymetrix platforms, GPL96 and GPL570, we defined stable gene pairs commonly detected by the two platforms as the ultimate stable gene pairs. Because different platforms have different probe designs and experimental protocol, some gene pairs may not keep consistent REOs in different platforms. Further study is needed to evaluate whether the two coupled signatures indentified in this study are suitable for microarray data produced by other platforms.

Data and pre-processing
All gene expression datasets for normal breast tissue and ER+ breast cancer were collected from GEO [61], as described in detail in Table 1. All samples used in this study fell into three categories: samples of normal breast tissue for identifying gene pairs with stable REOs in normal breast tissue, samples of ER+ lymph-nodenegative breast cancer patients accepting surgery only for developing a drug-free prognostic signature and samples of post-operative tamoxifen-treated ER+ breast cancer patients for developing a therapy benefit predictive signature. The third category included both lymph-nodenegative and lymph-node-positive patients, while most of them are in early stage (Table 2). RFS served as the prognosis endpoint, representing both disease-free survival and distant metastasis-free survival [47].
All the above-mentioned data were produced by the GPL96 or GPL570 platform. For each of the datasets, raw intensity files (.CEL) were processed using the RMA algorithm for background adjustment and median polish summarization without quantile normalization [62]. With the custom CDF file, each probe set ID was mapped to Gene ID, and then probe sets that mapped to multiple Gene IDs or did not map to any Gene ID were removed. The expression measurements of all probe sets corresponding to the same Gene ID were averaged to obtain a single measurement (on the log2 scale). The raw mRNA expression data of the post-operative tamoxifen treated patients were processed with the RMA quantile normalization algorithm in order to select DEGs between the high-and low-risk patients predicted by the drug-free prognostic signature.
The annotation data of 1320 canonical pathways, covering 8428 unique genes, were downloaded from the C2 collection of MSigDB (Version 4.0, updated May 31, 2013) [63] for personalized pathway analysis.

Consistency evaluation of stable REOs detected by different platforms
We focused on analyzing the 12752 genes measured by both the GPL96 and GPL570 platforms. For a collection of normal breast samples measured by a particular platform, if gene A had a higher (or smaller) Oncotarget 44605 www.impactjournals.com/oncotarget expression level than gene B in more than 99% normal samples, then the gene pair (A,B) was defined as stable gene pair. Based on the overlapping stable gene pairs detected by both the GPL96 and GPL570 platforms, a consistency score was calculated as the percentage of stable gene pairs with identical REOs in both collections of normal samples. We evaluated whether the consistency score was higher than what expected by chance using the binomial distribution test as following: where 0.5 is the probability of observing a gene pair having the same REO in two collections of normal samples by chance, n denotes the number of overlapping stable gene pairs detected by the two platforms, and k denotes the number of stable gene pairs with identical REOs in the two collections of normal samples.

Survival analysis
The univariate Cox proportional-hazards model [64] was used to evaluate the correlation of disruption indexes of pathways with the RFS and to evaluate whether a gene pair's reversal REOs were significantly correlated with poor RFS. When identifying RFS relevant gene pairs, we characterized REO of intra-pathway gene pairs for each sample as a binary vector in which 0 represented the REO of the intra-pathway gene pair in a cancer sample in line with that in normal tissue while 1 represented reversal REO. Kaplan-Meier survival plots and logrank tests [65] were used to evaluate the differences in RFS of distinct groups. The Cox proportional-hazards model was also performed to calculate the hazard ratios (HRs) and their 95% confidence intervals (CIs). The independent prognostic value of a signature was assessed by multivariate Cox proportional-hazards model. To evaluate the predictive performance of a signature we adopted the concordance index (C-index), which is a measure of overall concordance between predicted risk scores and observed RFS [66]. C-index, ranging from 0.5 (indicating random chance) to 1 (indicating perfect discrimination), is one of the most appropriate index for studies focusing on long-term risk prediction [67]. The Benjamini-Hochberg multiple testing correction was used to estimate the false discovery rate (FDR) [68]. All statistical analyses were performed using the R software package version 3.0.1.

Algorithm for searching optimum signatures
For a set of gene pairs whose REOs were associated with poor RFS, a forward-stepwise selection algorithm was performed to search for a optimal subset of these gene pairs that resulted in the highest C-index. Starting with the intra-pathway gene pair with the largest C-index as the seed signature, candidate intra-pathway gene pairs were added to the signature one at a time until the addition of one gene pair did not improve predictive performance.