We included 202 eligible patients from our previously published cohort of 465 primary TNBC patients treated at Fudan University Shanghai Cancer Center (FUSCCTNBC) (23). Patients were included based on the following criteria: histologic diagnosis of lymph node-negative TNBC with RNA-sequencing data and follow-up information for recurrence and metastasis. The RNA-sequencing data are available in the Sequence Read Archive (RNA-seq: SRP157974). Patients with contralateral breast cancer, lymph node recurrence, and unknown sites of recurrence were excluded. Lymph node status was independently confirmed by two experienced pathologists. The date of diagnosis of metastasis was defined when metastasis was either confirmed by biopsy or clinically diagnosed. The follow-up of this cohort was completed on June 11, 2019. Distant metastasis-free survival (DMFS) was defined as the interval between diagnosis and the first distant metastasis (viscera/bone/brain). Patients without events were censored from the time point of the last follow-up.

The present study was reviewed and approved by the Ethics Committee of Fudan University Shanghai Cancer Center (Ethics number: 050432-4-1212B). The patients provided written informed consent to participate in this study.

To identify mRNAs of prognostic value, analysis for differentially expressed mRNAs between two groups was performed using the limma package (version 3.48.0) in R software. We also performed Gene Set Enrichment Analysis (GSEA) of differentially expressed genes between the two groups with or without distant metastasis using the RNA-sequencing data and GSEA software (GSEA_4.1.0) (24, 25).

The cohort was randomly divided into the training set (n=142) and validation set (n=60) at a ratio of 7 to 3 by the caret package (version 6.0-88) in R software. Pearson chi-square test or Fisher’s exact test was used to ensure that there was no significant difference and that no bias was introduced in clinicopathological characteristics between the two sets. Least absolute shrinkage and selection operator (LASSO) Cox regression analysis was performed to further filter the differentially expressed mRNAs. A multivariate Cox regression model was used to determine the coefficient of each factor. The risk score of each model was used to estimate the probability of distant metastasis. The genomic risk score was calculated from individual gene expression measurements as follows: Genomic risk score = (β _B3GALT5-AS1 × B3GALT5-AS1) + (β _DNER × DNER) + (β _CSN1S1 × CSN1S1) + (β _KIF5A × KIF5A) + (β _SIX3 × SIX3) + (β _NOTUM × NOTUM) + (β _CPS1 × CPS1). The clinical risk score was calculated as follows: Clinical risk score = β_Age × Age (years)+ β_{Tumor size} × Tumor size (cm). The combined risk score was calculated as follows: Combined risk score = β_{Gene score} × Genomic risk score + β_{Clinical score} × Clinical risk score.

Patients were stratified into high- and low-risk groups based on optimum cutoff risk scores determined by the “surv_cutpoint” function in the survminer package (version 0.4.9) in R software. Kaplan-Meier analyses and log-rank tests were performed to assess the differences in DMFS between the high- and low-risk groups. The time-dependent receiver operating characteristic (ROC) curve was used to measure the prognostic performance by comparing the area under the ROC curve (AUC) values.

Based on data availability and clinical evidence (9, 26, 27), a nomogram was constructed integrating the seven-gene risk score, age of the patients at surgery, and pathological tumor size. We measured the predictive accuracy of the nomogram via Harrell’s concordance index (C-index) in the training and validation sets. In addition, the predictive capacity of the nomogram was also evaluated using calibration curve and decision curve analysis (DCA).

Pearson’s chi-square test or Fisher’s exact test was used to compare the clinical and pathological characteristics between the training set and validation set. All statistical analyses were performed using the SPSS 22.0 (SPSS Inc.) or R software (version 4.1.0, www.r-project.com). A value of P < 0.05 was considered statistically significant.

The clinical and pathological characteristics of 202 patients and their primary tumors are summarized in Table 1 . Of 202 lymph node-negative TNBC patients, the median follow-up was 68.2 months (interquartile range, 57.6-80.6 months). Overall, 12 (5.9%) cases with distant metastasis were observed. Of the 12 patients, 4 (33.3%) patients had multisite metastasis, and 7 (58.3%) patients died due to breast cancer during follow-up. The median tumor size and age of the patients at surgery in this study cohort were 2.5 centimeters (range 0.8-12.0) and 53 years (range 25-82), respectively.

An overview of the study design is shown in Figure 1 . Using log₂(fold change) > 1 or < -1 and P < 0.05, we identified 71 differentially expressed mRNAs between the two groups with or without distant metastasis. We also performed Gene Set Enrichment Analysis of differentially expressed genes between the two groups with or without distant metastasis using the RNA-sequencing data. In patients with distant metastasis, 25 gene sets were significantly enriched at nominal P value < 0.05. The top ten gene sets enriched in 12 lymph node-negative TNBC patients with distant metastasis compared to 190 patients without distant metastasis were illustrated in Figure S1 . In patients with distant metastasis, 56 mRNAs were upregulated, whereas 15 mRNAs were downregulated ( Figure 2 ). We constructed a matrix integrating RNA-sequencing data of 71 differentially expressed mRNAs and clinicopathological data of all 202 patients. Next, patients were randomly classified into the training set (n = 142) and validation set (n = 60). There was no difference in all characteristics between the training and internal validation sets ( Table 1 ). Seven genes, including B3GALT5-AS1, DNER, CSN1S1, KIF5A, SIX3, NOTUM, and CPS1, were selected using the LASSO Cox regression model in the training set. The summary of log₂(fold change), multivariable Cox regression coefficient, hazard ratio, 95% confidence interval, and P value for selected genes are presented in Table 2 . Time-dependent ROCs and Kaplan–Meier curves were used to evaluate the prognostic potential of the seven-gene signature for DMFS ( Figures 3A, B ). The AUC values for 3-, 4-, and 5-year DMFS were 0.823, 0.879, and 0.870 in the training set and 0.727, 0.705, and 0.689 in the validation set, respectively ( Figure 3A ). The formula of genomic risk score is as follows: genomic risk score = 0.18801037 × DNER + 0.28358112 × CSN1S1 + 0.36011127 × KIF5A + 0.57677377 × SIX3 + 0.70105693 × NOTUM + 0.74508978 × CPS1 - 0.06761698 × B3GALT5-AS1. Patients were stratified into high- (n = 15) and low-risk groups (n = 127) by selecting the optimal cutoff value (1.78) in the training set ( Figures 3B, C ). Using the same cutoff value (1.78), the patients were also divided into high-risk (n = 8) and low-risk (n = 52) groups in the validation set ( Figures 3B, C ). The Kaplan-Meier analyses for DMFS as a function of the seven-gene signature showed highly significant differences between the high- and low-risk groups ( Figure 3B , P < 0.001 in the training set; P = 0.039 in the validation set).

We also created a clinical prognostic model using the following clinically significant predictors: age and tumor size. The summary of multivariable Cox regression coefficient, hazard ratio, 95% confidence interval, and P value for age and tumor size are presented in Table S1 . The formula of clinical risk score is as follows: clinical risk score = 0.21532 × Tumor size (cm) - 0.04466 × Age (years). The AUC values of the clinical model for 3-, 4-, and 5-year DMFS were 0.755, 0.699, and 0.693 in the training set and 0.574, 0.651, and 0.631 in the validation set, respectively ( Figure 4A ). The genomic risk score remained an independent prognostic factor in the multivariate Cox analysis after adjusting for patient age and tumor size in both the training set (hazard ratio = 2.64, 95% CI: 1.76-3.96, P < 0.001) and validation set (hazard ratio = 1.63, 95% CI: 1.07-2.49, P = 0.02). The combined risk score was derived from the genomic and clinical risk score as follows: combined risk score = 0.9702 × Genomic risk score + 1.0854 × Clinical risk score. After integrating the clinical model with the genomic risk score, the AUC values for 3-, 4-, and 5-year DMFS were 0.836, 0.888, and 0.882 in the training set, respectively ( Figure 4B ). The AUC values of the combined model remained high in the validation set with values of 0.801, 0.793, and 0.768 for 3-, 4-, and 5-year DMFS, respectively ( Figure 4B ). Patients were stratified into high- (n = 15 or 9) and low-risk groups (n = 127 or 51) in the training set or validation set ( Figure 4C ). The Kaplan-Meier analyses for DMFS as a function of the combined model showed a significant difference between the high- and low-risk groups in the training set ( Figure 4C , P < 0.001). Likewise, the trend was also observed in the validation set ( Figure 4C , P = 0.071).

We integrated the seven-gene signature with age and tumor size to construct a prognostic nomogram in the training set ( Figure 5A ). The C-index value for the combined models was 0.874 in the training set and 0.805 in the validation set. The 4- and 5-year time-dependent ROC curves for the seven-gene, clinical, and combined models are illustrated in Figure 5B . Both the seven-gene model and combined model showed better prognostic performance than the clinical model for predicting 4-year DMFS (P = 0.046 for the gene model; P = 0.029 for the combined model). The combined model showed significantly better prognostic performance than the clinical model for predicting 5-year DMFS (P = 0.038), and the seven-gene model also trended toward significance (P = 0.065). The calibration analysis of the 4-year DMFS prediction is shown in Figure 5C . The solid blue line has a closer fit to the dotted gray line, indicating great predictive accuracy of the nomogram. Decision curve analysis (DCA) revealed that compared to the clinical model, the seven-gene model and combined model were superior in predicting 4-year DMFS ( Figure 5D ).