The following key words were searched for in the Gene Expression Omnibus (GEO) dataset: ‘Breast cancer’, ‘tamoxifen’, ‘expression’ and ‘microarray’; and datasets with <100 ER+ tamoxifen-treated samples, or datasets without survival information, were then manually filtered out. Following this, four datasets, GSE17705, GSE22219, GSE26971 and GSE56884, were then retained for further analysis. Furthermore, samples that were not primary tumor tissue were also excluded during this process. Raw data were then downloaded in the CEL format from the GEO datasets. Following this, background correction and normalization with Robust Multiarray Averaging were carried out using the R package ‘affy’ function ‘rma’ (v1.54.0). Probe and gene names were matched according to the manufacturer-provided annotation file. Genes with more than one complementary probe were merged and the average values were retained as the expression levels for the corresponding genes.

Cox univariate regression was implemented in both GSE26971 and GSE17005 datasets via correlation of each individual gene's expression with the survival information in both datasets using the R package ‘survival’ Genes significantly correlated with distant metastasis-free survival time in both GSE26971 and GSE17005 datasets were retained for further analyses as candidate genes. Random forest variable hunting was applied for the selection of a reasonable combination of candidate genes using R package ‘RandomForestSRC’ (v1.9.0). The parameter used was: 100 repeats and 100 iterations. Following this, multivariate Cox regression analysis was carried out in order to develop the linear risk score model using the selected candidate genes, and coefficients were solved with the training dataset, GSE17005. In the validation datasets (GSE22219, GSE26971 and GSE58644), these coefficients were locked in order to calculate the risk score of samples in the other datasets.

All statistical analyses were performed using R software (v3.0.1; https://www.r-project.org) and R packages. Normalizations of affymetrix raw data were performed with R package ‘affy’ using the function ‘rms’. The survival analysis and cox probability hazard model construction were performed with R package ‘survival’. Random forest variable hunting was implemented with R package ‘RandomForestSRC’, and receiver operating characteristic (ROC) curves were generated with R package ‘pROC’ (11). The nomogram was plotted with the clinical data in the training dataset using R package ‘rms’.

The detailed workflow of gene selection and model development is presented in Fig. 1A. The levels of association between gene expression levels and treatment outcomes (survival data) were assessed using Cox univariate regression. Genes associated with overall survival in both the GSE17705 and GSE26971 datasets were identified, and a total of 48 genes were then selected as candidates. Following this, the random forest variable hunting was performed in order to select for the optimal candidate genes. Following identification of 10 candidate genes (Fig. 1B), risk scores using Cox multivariate regression and expression of 10 genes were then calculated. The coefficients are presented in Fig. 1C, and parameters of Cox regression are shown in Table I. The risk scores were calculated as follows (where gene names represent their respective expression levels): Risk score = (0.299988203)*cyclin B2 (CCNB2) + (0.640775607) *glutamyl-prolyl-tRNA synthetase (EPRS) + (−0.756716676) *α-ketoglutarate dependent dioxygenase (FTO) + (0.117814961) *stem-loop binding protein (SLBP) + (0.245606283)*CTD small phosphatase 1 (CTDSP1) + (−0.161767842) *cyclin A2 (CCNA2) + (0.196307548) *bridging integrator 3 (BIN3) + (−0.618268545) *RNA binding protein with multiple splicing (RBPMS) + (0.580014194) *forkhead box D1 (FOXD1) + (−0.288974361) *gene encoding WD repeat of SOCS box containing 2 (WSB2).

Patients were divided into two groups, a high-risk group or a low-risk group, according to their median risk score. Following this, the difference in survival between the high-risk and the low-risk groups was calculated, and the results revealed that the high-risk group had a reduced relapse-free time compared with the low-risk group, with median survival times of 3.75 vs. 6.5 years, respectively (P<0.001; Fig. 2A). The high-risk group tended to represent early metastasis, and genes with high expression levels tended to have positive coefficients and genes with low expression tended to have negative coefficients (Fig. 2B). The 5-year distant relapse-free survival rate of the high-risk group was 75%; whereas this value was revealed as being 96% in the low-risk group. These results indicated that the developed risk score was an effective predictive indicator for the distant relapse survival time period of patients with ER+ breast cancer treated with tamoxifen.

Considering that the risk score staging system was developed based upon gene expression data in the GSE17705 dataset, there was a potential risk that the model would over-fit to the dataset. In order to assess the robustness of the developed risk score model, three independent datasets (GSE22219, GSE26971 and GSE58644) were used for further validation. Following the locking of the coefficients for the 10 genes, a risk score for each patient was calculated. In addition to patients belonging to the training dataset, the patients belonging to each of the three independent datasets were artificially divided into high-risk and low-risk groups using median risk score values as cutoff values. The patients with high-risk scores tended to have early relapse, as was similarly demonstrated in patients belonging to the training dataset (Fig. 3A). Furthermore, the gene expression profiles for the 10 genes in the both the low-risk and the high-risk groups resemble those demonstrated by the training dataset (Fig. 3B). These results demonstrate that the risk score model is robust across datasets for the prediction of distant relapse in patients with ER+ breast cancer treated with tamoxifen.

Subsequently, the associations between clinical parameters (stage, age, grade, lymph node invasion and primary tumor size) with the risk score were evaluated. As revealed in Fig. 4A, age, tumor stage and grade were significantly associated with the risk score; whereas the other clinical parameters were not (P>0.05). To facilitate the utilization of the risk score, a 5-year distant relapse nomogram was plotted (Fig. 4B). According to this nomogram, risk score was one of the most important metastatic indicators.