Prediction of Lymph Node Metastasis in Breast Cancer by Gene Expression and Clinicopathological Models: Development and Validation within a Population-Based Cohort Free

The benefit of extensive surgical excisions of axillary lymph nodes in primary breast cancer is being questioned on the basis of the advances in adjuvant therapy and the results of randomized trials on axillary management in sentinel lymph node–positive disease (1–5). Although these trials mark a growing interest in minimizing axillary surgery, sentinel lymph node biopsy (SLNB) is still performed in all clinically node-negative patients despite the notion that the majority of them have disease-free axilla (6, 7). For these patients, the invasive procedure has no therapeutic benefit and better tools to predict the tumor's metastatic capacity are needed for correct guidance on surgical and systemic treatments.

Five intrinsic molecular subtypes of breast cancer have been identified by gene expression profiles, for distinguishing phenotypes and prognosis (8–10) and for guiding decisions on adjuvant treatment (11, 12). Although the prognostic value of gene profiling has been validated in node-negative disease (13), its prognostic utility has been extended to node-positive cases (14, 15). Despite the accumulated evidence of genomic profiles as prognostic indicators, nodal status retains its importance as an independent prognostic marker (16). Accurate preoperative assessment of axillary status and improved understanding of the biological processes associated with nodal metastasis are thus important for guidance on breast cancer treatment. Histopathological features, such as tumor size and lymphovascular invasion (LVI), have proved essential for assessing risk of nodal involvement and have been integrated into nomograms and risk models, together with multifocality and age among other features (7, 17). The first nomogram was developed in 2007 at the Memorial-Sloan Kettering Cancer Center (MSKCC) with an AUC of 0.75 and has been validated in independent cohorts, including European and Chinese populations (18–20), with AUCs ranging from 0.67 to 0.78. However, risk assessment of metastatic axillary spread based on these clinicopathological variables is considered imperfect and consequently SLNB remains the standard axillary staging procedure (18, 21, 22). Previous studies of gene expression patterns from the primary breast tumor for prediction of nodal metastasis have shown inconsistent results, ranging from being nearly-perfect to almost inadequate in small, selected cohorts (23–28). The classifiers developed to date are not sufficient to support relevant clinical conclusions (29). A remaining question is thus whether gene expression profiling alone or in combination with clinical variables may surpass the performance of existing methods, for instance, nomogram methods that report AUCs between 0.70 and 0.75, in performance (18, 30).

In the current study, we aimed to evaluate the potential of clinicopathologic and gene expression-based predictors (single and mixed) of lymph node metastasis in a population-based contemporary breast cancer cohort (n = 3,023) analyzed by RNA sequencing (RNA-seq) to identify patients least likely to benefit from surgical axillary staging, and consequently reduce the rate of unbeneficial SLNB. Furthermore, we aimed to evaluate biological processes related to nodal metastasis and estimate the prognostic impact of these predictors.

The study was approved by the Regional Ethical Review Board in Lund, Sweden (Registration numbers 2009/658, 2009/659, 2010/383, 2012/58, 2012/379, 2013/12, 2013/459, 2014/521, 2014/681, 2015/277, 2016/541, 2016/944, 2018/267) and was conducted in accordance with the Declaration of Helsinki. All patients were included after obtaining written informed consent.

The Swedish National Quality Registry for Breast Cancer (31) identified 5,235 patients undergoing primary surgery in the Southern Healthcare Region from September 1, 2010, to March 31, 2015. Of these, 4,633 patients were enrolled in the Sweden Cancerome Analysis Network–Breast (SCAN-B). SCAN-B employs a prospective, population-based, observational cohort study design (ClinicalTrials.gov ID: NCT02306096; refs. 32, 33) and the option for inclusion was proposed to all patients with primary breast cancer in the catchment area. High-quality tumor RNA-seq profiles were obtained for 3,023 patients (Fig. 1A). Information on age, tumor size, nodal status, mode of detection, multifocality, histologic grade, LVI, status of estrogen receptor (ER), progesterone receptor (PR), and HER2, adjuvant therapies, and overall survival (OS) were retrieved from the Swedish National Quality Registry for Breast Cancer (31) and Statistics Sweden (34). Ki-67 values were excluded due to incomplete data. Exclusion criteria for patients included neoadjuvant therapy and unknown or inconsistent pathologic nodal status from cancer registry data (Fig. 1A). Tumors were classified into four surrogate molecular subtypes based on the ER, PR, and HER2 status registered in the Swedish National Quality Registry for Breast Cancer: (i) all cases, (ii) ER⁺HER2⁻, (iii) HER2⁺, and (iv) triple-negative breast cancer (TNBC). Patients with micro- or macrometastases on SLNB were defined as lymph node-positive and were recommended completion axillary lymph node dissection (ALND) according to the Swedish National Guidelines for Breast Cancer at that time (35). The median follow-up for patients without an event was 4.53 years (0.16–6.78).

Gene expression profiling was performed using RNA-seq as previously described (32, 33) and reported elsewhere (GSE96058; ref. 36) with complementary material in Supplementary Table S1. In brief, fresh breast tumor specimens, taken by a breast pathologist, were stored cooled in microtubes with the preservative solution (RNAlater) before delivery to the central SCAN-B laboratory in Lund, Sweden for tissue sample treatment and RNA-seq. Sequencing was performed to a depth of approximately 30 million paired-end reads (Illumina HiSeq2000 or NextSeq500 systems) using a dUTP protocol. The quality of RNA-seq data were assessed; aligned reads < 10 million, duplication level > 55%, or RNA integrity number (RIN) <6 were excluded. A gene was represented by a single combined FPKM value of the matching transcript. Only refseq annotated genes (n = 19,102) were used for further analysis. An offset of +1 was added to all FPKM values followed by log₂ transformation. Intrinsic breast cancer subtype classification was performed by applying absolute intrinsic molecular subtyping (AIMS; ref. 37) to assign the PAM50 subtypes. AIMS subtype classification can be accomplished by a gene expression profile for each sample without the requirement of a normalized dataset. Further details on performed calculation of biological metagene scores (38) can be found in Supplementary Information, Appendix S1. Pathway analysis of differentially expressed genes was performed using the topGO R package as described (Supplementary Information, Appendix S1).

The 3,023 patients were divided into development (n = 2278) and validation (n = 745) cohorts. The validation cohort comprised of patients diagnosed in 2011 (>5-year follow-up), while the remaining cases were assigned to the development cohort. In the development cohort, seven machine-learning models were applied to each of the four evaluation groups. For each evaluation group and model type, three classifiers were trained based on (i) clinicopathologic features alone (CLINICAL), (ii) RNA-seq data alone (GEX), and (iii) mix of clinicopathologic features and RNA-seq data (MIXED), totaling 4 × 7 × 3 models. To identify the machine learning model with the best performance, the evaluation groups of the development cohort were divided into internal training and test sets by splitting the data into either 80% training and 20% test, or 50% training and 50% test (Fig. 1B). The training:test ratios were chosen based on the differences in sizes of the evaluation groups, which was due to the population-based nature of the SCAN-B cohort. Internal validation was performed by 10 iterations of 4-fold cross-validations, repeated 10 times (Fig. 1B). Assessment of model performance (the area under the ROC curve, AUC) identified generalized boosted regression models (GBM) with the highest predictive performance across evaluation groups and model types (CLINICAL, GEX, MIXED; Supplementary Fig. S1). A final GBM model (CLINICAL, GEX, MIXED) was then retrained using all samples in the development cohort for each evaluation group, removing features with no estimated importance (totaling 4 × 3 final GBM models). The inclusion of all samples in the development cohort for retraining the final locked models was to provide as many training cases as possible for GBM, as well as to avoid selectively choosing one model from the iterative internal validation process. Genes were eligible for training the final GEX and MIXED models if they were among the 5,000 most varying for each evaluation group and had a t test derived P value for difference between N0 and N+ samples <10% after Benjamin–Hochberg correction for multiple testing. Applying this approach, 74 to 1,633 eligible genes, depending on classifier and evaluation group, were included in the final training procedure (Supplementary Fig. S2). The final models were applied to the independent validation cohort for each evaluation group. Methodological details of the training procedure, feature selection, and deployed models are provided in the Supplementary Information, Appendix S1. Final data object models, including prediction scripts, are available on request from the authors.

A cutoff to discriminate nodal-negativity was set to correspond to maximized negative predictive value (NPV). This level of cutoff was intended to identify individuals with a very low probability of axillary disease where surgical axillary staging by SLNB would not be supported by the models. True positive (TP), true negative (TN), false positive (FP), and false negative (FN) results were assessed. The false negative rate (FNR) was calculated as the number of false negative cases divided by the number of all cases with axillary nodal metastasis (FNR = FN/FN+TP). This portrays the paucity of the models in predicting nodal spread as a ratio relative to all pathology-verified axillary nodal disease burden. Possible SLNB reduction rate = (TN+FN)/(TN+FN+TP+FP) was assessed at a cut-off equivalent maximized NPV and at alternative cut-offs corresponding to the maximum FNRs of 5% and 10% reflecting accepted FNR of the SLNB procedure.

Survival curves were estimated by the Kaplan–Meier method and compared using the log-rank test. OS was the endpoint for survival analyses. HRs were estimated by univariable or multivariable Cox regression. SPSS Statistics for Windows version 24 (IBM Corp.) was used for descriptive statistics with Pearson χ² test. Machine learning and survival analyses were performed in R (39) using the Caret package and the coxph R function. The developing of the predictive models and the reporting of the findings were in accordance with an EQUATOR Guideline for reporting machine learning predictive models (40) and STARD 2015 checklist (41).

The clinicopathologic characteristics for all patients in the catchment region and among those enrolled in SCAN-B, both in the development (n = 2,278) and validation cohorts (n = 745) were comparable (Table 1; Supplementary Table S1, page 2). However, tissue for RNA-seq was not always obtainable for tumors ≤10 mm. SLNB alone, without completion of ALND, was performed in 62% of patients. Axillary metastasis, N+, was found in 36% and 34% of the patients in the development and validation cohorts, respectively. Information on the metastatic size was incomplete; micrometastasis was the largest deposit in 6% of node-positive patients with obtainable records. There were significant differences in the proportion of lymph node metastasis across AIMS-PAM50 molecular subtypes within the development (P <0.001) and validation cohorts (P = 0.009; Supplementary Table S2).

Seven different machine-learning methods were evaluated and GBM was selected on the basis of best performance (Fig. 1B). GBM, in the implementation used by us, closely follows Friedman gradient boosting machine (42) and involves a stochastic gradient boosting strategy for logistic regression using a Bernoulli distribution (due to binary outcome of N0/N+). Gradient boosting combines multiple simple models “weak learners” into a single “strong” composite model through an iterative process. Figure 2A summarizes performances of the GBM models for prediction of nodal metastasis in the complete development cohort, with MIXED predictors achieving highest performance in all evaluation groups. Figure 2B shows predictive performance in the validation cohort. The overall levels of accuracy in the independent validation cohort suggest model robustness. However, a slight drop in performance was observed for HER2⁺ and TNBC validation groups for predictors including gene expression features, which could indicate possible overfitting in the training procedure as well as a reflection of small subgroups. Although the best predictive performance for the total validation cohort (n = 745) was achieved by a MIXED predictor [AUC = 0.72, 95% confidence interval (CI) = 0.68–0.76], the CLINICAL predictor showed comparable performance (AUC = 0.71; CI = 0.67–0.75). For the ER⁺HER2⁻ validation group, the CLINICAL model achieved an AUC of 0.73 (CI = 0.68–0.77), while the MIXED predictor achieved an AUC = 0.72 (CI = 0.68–0.77). For the HER2⁺ validation group, the CLINICAL predictor showed an AUC = 0.70 (CI = 0.58–0.81), and the MIXED model showed an AUC = 0.66 (CI = 0.55–0.78). In the TNBC validation group, the CLINICAL and MIXED predictors achieved AUCs of 0.65 (CI = 0.53–0.76) and 0.68 (CI = 0.56–0.81), respectively. Discriminatory performances of the GEX predictors were inferior to those of the CLINICAL and MIXED models, with AUCs ranging from 0.58 to 0.67. Although discriminatory performances were highly similar between the CLINICAL and MIXED predictors in general, Cohen kappa coefficients were 0.54, 0.56, 0.61, and 0.31 for the overall validation cohort, ER⁺HER2⁻, HER2⁺, and the TNBC subgroups, respectively (Supplementary Fig. S3A). These values indicate a fair to moderate agreement between the two models. In comparison, better agreement was observed between the GEX and MIXED models in TNBC tumors (kappa coefficient = 0.66), and a modest agreement for the remaining groups (Supplementary Fig. S3B). When comparing the agreement between predictors (i.e., proportion of cases predicted similarly by two predictors), we found that 75% to 85% of cases were predicted similarly for the two different predictor set comparisons across evaluation groups (Supplementary Fig. S3C). Finally, when stratifying N+ disease (the training label) to N1 (1–3 metastatic nodes) and N2 (≥4 metastatic nodes), we noted that most discrepant classification between CLINICAL, MIXED, and GEX predictors and the pathologic class occurred in the N1 group across all evaluation groups (Supplementary Fig. S3D). This observation indicates that N1 tumors have a heterogeneous phenotype in terms of the investigated clinicopathologic variables and gene expression patterns.

Table 2 displays possible SLNB reduction rates after applying the prediction models to assign disease-free axilla (N0 vs. N+) in the SCAN-B validation cohort (both for the complete SCAN-B 2011 validation dataset and the largest evaluation group comprising ER⁺HER2⁻ tumors within the validation dataset) corresponding to cutoffs at the maximum negative predictive value (FNR 0.4%) and FNR 5% and 10% (Fig. 2C). In the overall validation cohort, a SLNB reduction rate of 2.8% would be attained by applying CLINICAL predictors at a cut-off corresponding to maximized NPV 95% (FNR 0.40%) to discriminate disease-free axilla. Allowing a maximum of 5% to 10% FNR in discriminating N0 by CLINICAL predictors, possible SLNB reduction rates of 11% to 23% were displayed. However, if MIXED predictors were applied, SLNB reduction rates of 19% to 30% would be yield at a maximum of 5% to 10% FNR (corresponding to the accepted level of FNR for the SLNB procedure). In general, MIXED predictors yielded 6% to 7% further reduction rates compared with CLINICAL models at 5% to 10% FNR.

Variables in the final models are listed in Supplementary Table S3 for each predictor type and evaluation group. Among the CLINICAL predictors, the most significant variables associated with axillary metastasis were tumor size, LVI, age, and multifocality. Tumor size and LVI remained the top two variables associated with nodal disease for MIXED predictors in the validation cohort for all cases, ER⁺HER2⁻ cases, and HER2⁺ cases, but ranked seventh and sixth for the TNBC group, respectively.

Biological processes related to predictor classification were assessed by metagene expression analysis in four nodal status subgroups defined by combinations of pathologic anatomic diagnosis of the nodal status, N0/N+ (PAD−/PAD+), and model predicted axillary status, N0/N+ (Pred−/Pred+; Fig. 3A–C). For the ER⁺HER2⁻ cases, differences in expression of proliferation-related genes (checkpoint metagene) were noted among Pred−/Pred+ samples, regardless of the pathologic status (Fig. 3A).

PAM50 classification revealed heterogeneity in the four nodal status subgroups in the ER⁺HER2⁻ evaluation group (Fig. 3D). Nevertheless, distribution of the PAM50 subtypes was comparable in a pair-wise pattern defined by predictor nodal status (Pred−/Pred+) rather than by pathologic status (PAD−/PAD+) and involved higher proportions of the luminal B intrinsic subtype in the predicted N+ cases. Similar patterns were observed in the HER2⁺ evaluation group; samples with model-predicted nodal metastasis showed higher expression of proliferation-related genes (Fig. 3B) and intrinsic molecular subtypes associated with a more aggressive tumor type (Fig. 3E). For the TNBC evaluation group, low expression of basal-like markers was detected in samples with predicted nodal metastasis by the GEX and MIXED models but not by the CLINICAL model (Fig. 3C).

Unsupervised hierarchical analysis of the validation cohort was performed to investigate the relationship between global gene expression patterns and the distribution of nodal status, clinicopathologic variables, and breast cancer subgroups. Fig. 3G shows the transcriptional landscape in the context of N0/N+ status. Heterogeneous distribution was observed for N0/N+ status in relation to key clinicopathologic variables and breast cancer subgroups, irrespective of whether these subgroups were defined by surrogate subtype markers or by molecular PAM50 classification. Functional analysis based on gene ontologies (GO) for derived GEX and MIXED predictors are available in Supplementary Fig. S2, demonstrating strong association with proliferation related biological processes especially in the ER⁺HER2⁻ group, but also in the HER2⁺ group.

Assessment of prognosis based on nodal status was performed in patients with ER⁺HER2⁻ disease (n = 398) who received adjuvant endocrine therapy only (the largest subgroup). Other treatment groups, such as chemotherapy plus endocrine therapy, did not include enough patients to allow statistical analysis and evaluation at this time. Figure 4A presents OS for the cohort stratified by the combinations of pathologic and predicted axillary status (PAD−Pred−, PAD+Pred−, PAD−Pred+, PAD+Pred+). On the basis of the available follow-up data, positive pathology-based lymph node status did not predict a significantly worse OS in Kaplan–Meier analysis (log-rank P = 0.13, Fig. 4A) or univariate Cox analysis (HR = 1.58; 95% CI = 0.87–2.86; Fig. 4B), although the HR indicates a trend toward significance. In contrast, predicted nodal status stratified patients into better (PAD−Pred− and PAD+Pred−) or worse (PAD−Pred+ and PAD+Pred+) outcomes irrespective of pathologic status, for all models. In univariable analysis, factors predicting OS included age (years), tumor size (mm), histologic grade 3 versus grade 1 and grade 2, node positivity versus node negativity predicted by any of the models, and combinations of the pathologic and predicted axillary status (PAD+Pred− or PAD−Pred+ or PAD+Pred+ versus PAD−Pred−) across models. Next, we tested whether N0/N+ status by PAD added prognostic information in multivariable Cox regression models with the predictions from GEX, CLINICAL, or MIXED as the only covariates included in the model, respectively. However, in none of the three analyses N0/N+ by PAD reached statistical significance (P > 0.05), while the GEX, CLINICAL, and MIXED predictions did (P < 0.05). It should be noted that our derived models were not optimized for outcome; still they capture a signal showing prognostic value irrespective of the pathologically defined nodal status. Finally, the prognostic importance of the PAD+Pred+ versus the PAD−Pred− remained significant (HR = 2.13; CI = 1.06–4.26; P = 0.03) in a Cox multivariable analysis based on the GEX predictor when adjusted for tumor size (mm) and age (years) (Fig. 4B and C), but was not significant (HR = 2.04; CI = 0.93–4.49; P = 0.07) when adjusting for histologic grade. CLINICAL and MIXED predictors were not tested in the multivariable analysis as they included the adjusted covariates.

This study presents data from the largest prospectively collected and population-based cohort to date on genomic predictors for nodal metastasis, and comprehensively evaluates and independently validates mRNA-seq and clinicopathologic features as predictors of nodal metastasis in ER⁺HER2⁻, HER2⁺, and TNBC tumors. Today (43), a precise noninvasive option for nodal staging would circumvent the adverse effects of axillary surgical staging since the majority of breast cancer patients present with low-risk tumors, node-negative disease, and excellent overall prognosis. In the overall validation cohort, the predictor using RNA-seq features and clinicopathological characteristics (MIXED) showed the highest accuracy in predicting nodal metastasis, attaining AUC = 0.72. Thus, the ability to classify nodal status using current predictors and endpoint definition was not perfect. Nevertheless, the nodal metastatic features recognized by these predictors implied an added prognostic impact compared to pathological nodal status alone and the presented SLNB reduction rate was higher for the MIXED predictor compared to the CLINICAL predictor. Notably, the AUCs attained in our models are similar to those of prediction models previously published for nodal metastasis (AUC 0.67–0.78), which included only clinicopathologic data (18–21, 30), but did not report on possible SLNB reduction rate.

Although the addition of gene expression data to existing clinicopathologic variables did not show a clear superiority in predicting nodal status, it should be noted that machine-learning approaches identify the most dominant trait for each predictor type and capture different metastatic features. Even though the discriminatory performances (AUC) of the CLINICAL and MIXED predictors were identical, there was only a moderate agreement between them in the validation cohort (Supplementary Fig. S3). Previous studies investigating prediction of nodal metastasis in breast cancer have reported diverse performances of GEX-based predictors, with AUCs between chance level to near perfect separation (23–28), likely to be due to differences in patient characteristics, cohort sizes, definition of nodal disease, gene expression analysis platforms, and feature-selection strategies. In other malignancies, nodal prediction accuracies between 50% to 60% have been reported based on TCGA data (28). Our predictors containing gene expression data (GEX and MIXED models) revealed AUCs of 0.58 to 0.72 during validation and were consistent with a mid-range performance as defined in previous publications (23–28). The current results are, however, validated in an unselected contemporary large primary breast cancer cohort. Still, although the small proportion of HER2⁺ and TNBC tumors reflect the breast cancer population in Sweden (31), results from these evaluation groups must be interpreted with caution.

One problem with developing mixed gene expression classifiers for the classification question at hand is the scarcity of larger cohorts including gene expression- and relevant clinicopathologic data that is commonly used in, for example, nomogram methods (such as lymphovascular invasion, mode of detection, multifocality and tumor location; refs. 18, 30). To further validate the MIXED classifier for ER⁺HER2⁻ tumors, we used the TCGA dataset (44) comprising 503 primary breast tumors. Despite considerable imputation of missing values for the important clinicopathologic parameters and caution must be exercised when interpreting these results, we found an AUC of 0.68 compared with 0.72 in our SCAN-B 2011 validation cohort. This is in line with the decline in performance reported for the MSKCC nomogram in an independent cohort (Supplementary information, Appendix S1; ref. 20). Corresponding possible SLNB reduction rates for patients with ER⁺HER2⁻ tumors within the TCGA cohort would be 12% to 18% if applying the MIXED predictor and accepting false negative rates of 5% to 10%. Based only on the AUC values for N0/N+ prediction by MIXED and CLINICAL models, differences appear too small to justify the added use of RNAseq. However, when considering the SLNB reduction rate as well, there is a clear indication that the MIXED model is able to identify more node negative patients than the CLINICAL model (exemplified by 7% absolute improvement in reduction rate for patients with ER⁺HER2⁻ tumors while accepting a false negative rate of 5%). With declining prices of GEX assays and the potential to combine different prognostic or predictive gene signatures in one transcriptional analysis, even a minor improvement of SLNB reduction rate in favor of the MIXED predictor can potentially be translated into cost effectiveness in the future.

Although the discrimination of nodal metastasis is challenging due to the heterogeneity of breast cancer, AUC analysis allows for different cutoffs to be assigned, determined by the prevalent clinical setting (45). For the purpose of identifying patients with very low risk of nodal metastasis to omit surgical axillary staging (SLNB), cutoffs were set to achieve the highest possible negative predictive value (Max NPV) in the prediction of N0. Moreover, additional cut-offs were set to reflect accepted FNR of the SLNB procedure (FNR 5%–10%) in itself, which is the golden standard for axillary staging. Allowing FNR 5% to 10% for N0 prediction, CLINICAL predictors and MIXED predictors would yield SLNB reduction rates of 11% to 23% (FNR 5%) and 19% to 30% (FNR 10%), respectively, supporting that the MIXED predictors identified more node-negative patients compared with the CLINICAL predictors. The importance of gene expression profiles for adjuvant treatment recommendation is today readily accepted, but not included as a tool for selection of optimal surgical strategy (11, 12, 15). If prospectively validated, integration of gene signatures in the preoperative diagnostic workup might reduce the rates of SLNB for patients at low risk of nodal involvement.

In line with reports on clinicopathologic variables associated with nodal metastasis (18, 21, 46, 47), tumor size and LVI were key variables in our CLINICAL and MIXED models for the ER⁺HER2⁻ and HER2⁺ groups, but not for TNBC tumors. Interestingly, TNBC is characterized by a low risk of nodal metastasis and displays a non-linear association between tumor size and the risk of nodal engagement (48, 49). We observed low expression of basal-like markers for TNBC tumors with predicted nodal metastasis by the GEX and MIXED models, but not for the CLINICAL model. In accordance with TNBC showing the lowest prediction overlap between CLINICAL and MIXED models (Supplementary Fig. S3A), five genes ranked higher in the MIXED TNBC model (MAGEA4, FUT3, MIA, LPAR2, LOXL4) than did LVI and tumor size. Interestingly, LOXL4 expression has been shown to be associated with tumor growth and lung metastasis in TNBC (50). The fair to moderate overlap between predictions by the CLINICAL, MIXED, and GEX models derived from the same samples indicates that these models capture different aspects of the N0/N+ disease, exaggerated in TNBC tumors when including gene expression patterns. Consequently, a great molecular heterogeneity must exist in each evaluation group given the chosen nodal classification endpoint.

To test the hypothesis of heterogeneity within the N0/N+ endpoint, we performed an unsupervised analysis of the validation cohort, demonstrating transcriptional heterogeneity of axillary status across the molecular breast cancer subtypes (Fig. 3G). Specifically, the N0/N+ status may not significantly reflect transcriptional variations in breast cancer and consequently presents a challenge for expression-based predictors. To substantiate this claim, we used a machine-learning procedure that included a feature selection step to identify genes differentially expressed between N0 and N+ in each evaluation group, with adjustment for multiple testing. For the MIXED classifiers this analysis identified 1268 differentially expressed genes in ER+HER2- tumors, 74 genes in HER2⁺, and 130 genes in TNBC (Supplementary Fig. S2). The low numbers of differentially expressed genes in the HER2⁺ and TNBC evaluation groups did not support a strong gene expression signal. Additionally, functional pathway analysis of the predictor gene sets showed heterogeneity in enriched terms, consistent with the observed transcriptional heterogeneity of the N0/N+ status. This finding is not surprising due to the limited number of genes retained after the retraining step (Fig. 1B) and the feature selection/optimization process of machine learning methods. However, for the GEX and MIXED predictor gene lists prior to the retraining step that excluded genes with weight ≤0 for the final locked models, we found a notable enrichment of proliferation-associated gene ontology terms in the ER⁺HER2⁻ and HER2⁺ subgroups (Supplementary Fig. S2). This is consistent with (i) the differential expression of the proliferation metagene, stratifying tumors based on predictor classifications rather than on pathologic status, and (ii) the variations in proportion of luminal B classified tumors for the same groups (Fig. 3), which could explain the observed prognostic disparities (Fig. 4).

The recognition of proliferation-related prognostic transcriptional programs, (51–55) especially in luminal breast cancer, prompted us to investigate the predictors' prognostic performance (for which they were not trained) in the largest subgroup in the validation cohort: ER⁺HER2⁻ tumors treated with adjuvant endocrine therapy only. In line with the expression pattern of the proliferation metagene (Fig. 3), the poorest outcome was observed for patients with predicted nodal disease involvement (Pred+), with worst OS for those with both pathologic and predicted nodal positivity (PAD+Pred+; Fig. 4A and B). This suggests that at least for patients with ER⁺HER2⁻ tumors PAD does not add prognostic value when nodes are predicted by the models to be negative. The weak association between pathology-based lymph node status and outcome is likely due to the short follow-up time and use of OS as clinical endpoint. In multivariable analysis adjusted for tumor size and age, the GEX PAD+Pred+ group remained significant. However, this did not prevail after further adjustment for grade. Therefore, the prognostic impact obtained from the predictors' proliferative features is related to the prognostic importance of histological grading. The limited cases within other evaluation groups precluded evaluation.

The current study has a number of limitations. First, clinicopathologic characteristics were obtained from a national registry, with incomplete data for some variables that could otherwise be included in predictor training or be used in defining evaluation groups, for xample, surrogate molecular subtype or Ki67 scores. Moreover, despite analyzing >3,000 cases, the number of HER2⁺ and TNBC tumors in our population-based cohort was limited and precluded assessment in more defined subgroups (e.g., HER2⁺ER⁻). Consequently, it may be hypothesized that a more refined molecular subgrouping of breast cancer by means of intrinsic signatures could improve the performance of gene expression-based classifiers. Second, our outcome analyses are limited by the type of endpoint available with complete coverage (OS) and the relatively short follow-up time of the prospective SCAN-B cohort. Finally, although we exhaustively evaluated seven different machine-learning models and different feature-selection approaches, we acknowledge that additional optimization and combinations of clinicopathological and transcriptional data could improve the results. Although the latter may be technically achievable, a highly relevant issue in moving gene expression predictors towards the clinic is the ability to predict single samples preferably without the need of reference sets, representing a limitation of our current models. Here, reports of methods deriving gene rule-based predictors (37, 56) may represent an alternative machine-learning approach, provided sufficient performance and adequate training numbers (57).

In summary, this study provides a first estimate of the performance of predictors including RNA-seq data as well as clinicopathologic data of axillary nodal status in a population-based setting and proposes possible SLNB reduction rates by applying predictors to discriminate N0. CLINICAL and MIXED predictors of nodal metastasis had comparable accuracy in line with previously published predictors, while more node-negative patients were identified by the MIXED predictors. The GEX models had inferior predictive performance, underscoring that the complexity of lymphatic spread dissemination is not only related to gene expression of the primary tumor. To fully harness the power of transcriptional profiling for nodal prediction, a more precise stratification of breast cancer subgroups also appears relevant. Given such improvements, predictors such as the MIXED classifier, combining clinicopathologic variables and gene expression data may develop into clinically relevant lymph node predictor with putative prognostic value.

A. Ehinger reports receiving speakers bureau honoraria from Roche, Amgen, and Novartis. Å. Borg reports receiving speakers bureau honoraria from AstraZeneca and Roche. No potential conflicts of interest were disclosed by the other authors.

Conception and design: L. Dihge, P.-O. Bendahl, Å. Borg, J. Staaf, L. Rydén

Development of methodology: L. Dihge, J. Vallon-Christersson, J. Staaf, L. Rydén

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): J. Vallon-Christersson, C. Hegardt, L.H. Saal, C. Larsson, A. Ehinger, N. Loman, M. Malmberg, L. Rydén

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): L. Dihge, J. Vallon-Christersson, L.H. Saal, J. Häkkinen, M. Malmberg, P.-O. Bendahl, J. Staaf, L. Rydén

Writing, review, and/or revision of the manuscript: L. Dihge, J. Vallon-Christersson, C. Hegardt, L.H. Saal, J. Häkkinen, C. Larsson, A. Ehinger, N. Loman, M. Malmberg, P.-O. Bendahl, Å. Borg, J. Staaf, L. Rydén

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): J. Vallon-Christersson, C. Hegardt, J. Häkkinen, M. Malmberg, J. Staaf, L. Rydén

Study supervision: J. Staaf, L. Rydén

The authors would like to acknowledge patients and clinicians participating in the SCAN-B study, all personnel at the central SCAN-B laboratory at the Division of Oncology and Pathology, Lund University, the South Swedish Breast Cancer Group (SSBCG), the Swedish National Breast Cancer Quality Registry (NKBC), The Regional Biobank Centre (RBC) Syd and the Regional Cancer Centre South (RCC Syd). This study was supported by grants from The Swedish Cancer Society (CAN 2016/659, CAN 2018/685), The Skåne County Councils Research and Development Foundation (2016/623611), The Governmental Funding of Clinical Research within the National Health Service (ALF) (2014/434901, 2018/40612, 2019/40304), The Swedish Breast Cancer Association (2017/1213), Mrs. Berta Kamprad Foundation (FBKS 2018-3-166), Mats Paulsson Foundation for Research, Innovation and Community Development (IACD 2017), The Crafoord Foundation (20180543), and The Erling-Persson Family Foundation (2018/1008). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.