Three independent patient cohorts with available RS were collected from two institutions and divided into training and validation sets based on the data sources and WSI resolutions. RS was defined as low (≤15), intermediate (16–25), and high (26–100) according to the results from the TAILORx trial (30). ER, PR, and HER2 interpretations were based on the updated ASCO/CAP recommendations (31, 32). All patients received surgery.

Training set: A total of 125 cases of ER+/HER2-/LN- breast cancer with RS diagnosed from 2011 to 2015 were collected from the Ohio State University. The RS ranged from 0 to 40. Among these 125 cases, 53, 59, and 13 cases had low scores, intermediate scores, and high scores, respectively.

Validation set 1: A total of 82 cases of ER+/HER2-/LN- breast cancer with RS diagnosed from 2012 to 2014 were retrieved from the Emory University. The RS ranged from 0 to 52. Among 82 cases, 40, 15, and 27 cases had low scores, intermediate scores, and high scores, respectively.

Validation set 2: Additional 175 cases of ER+/HER2-/LN- breast cancer with RS diagnosed from 2012 to 2014 were retrieved from the Emory University. The RS in this dataset ranged from 0 to 100. Among 175 cases, 68, 73, and 34 were low-, intermediate-, and high-score cases, respectively.

All three datasets included age at diagnosis, ER and PR IHC staining percentage (0–100) and intensity (1, 2, and 3), HER2 amplification by IHC and FISH (negative and equivocal), Nottingham tumor grade, and tumor size. Additional features retrieved for validation sets 1 and 2 included Ki-67 score, stage, chemotherapy, radiation therapy, overall survival (OS), disease-free survival (DFS), and distant metastasis (metastasis other than axillary LN metastasis). One representative tumor hematoxylin and eosin (H&E) stained WSI from each case in the training set and validation set 1 was scanned at 40 × magnification and validation set 2 at 20 × with an Aperio AT2 scanner.

The clinicopathological information of these three datasets is summarized in Table 1. The ER and PR expressions for all three cohort datasets were evaluated with an H-score (percentage × intensity). This study was approved by the Institutional Review Board at the Emory University and the Ohio State University.

Image normalization: As 40 × images have a higher resolution for annotations, we chose the 40 × for data analysis. After linearly resizing with a scaling factor of two along the image width and height directions, all images in validation set 2 had the same magnification of 40 × as training set and validation set 1. We also used the sparse non-negative matrix factorization-based color transfer method (33) to normalize the image color styles in all three cohort datasets (Figure 1).

Data preprocessing for DL training: Although we had three datasets for RS prediction analysis, we used two independent image datasets for cell detection and segmentation training, one from our lab and the other from the public MoNuSeg-2018 dataset. We collected 797 images with tumor nuclei point annotations, 500 images with TIL point annotations, and 26 images with annotations of tumor nuclei contours from the independent dataset. All the annotations were produced and confirmed by the pathologists (Supplementary Figure 1). Two pathologists made the annotations with Aperio ImageScope and GIMP. Additionally, 30 H&E images from the public MoNuSeg-2018 dataset were used in the segmentation dataset. They had annotations of cell nucleus contours (Supplementary Figure 1C). Each DL dataset was randomly divided into training, validation, and testing groups with an approximate proportion of 70:15:15.

For detection, classification, and segmentation analyses, we used the Mask R-CNN (MRCNN) (34) to construct the image processing models in this project. MRCNN was extended from Faster R-CNN (35) that was in turn developed based on Fast R-CNN (36). The overall schema of the developed WSI image processing pipeline is presented in Figure 2. The DL MRCNN pipeline was constructed with library TensorFlow and Keras. The image processing module contained three MRCNN models specifically for tumor cell detection, TIL detection, and tumor nucleus segmentation, respectively. Image tiles with tissue were extracted from WSIs by thresholding the “Saturation” channel of the HSV color space with the threshold set to 30. Each image tile was then analyzed by three MRCNN models separately. The center of each bounding box is considered the center of a detected cell of interest. The segmentation branch in the MRCNN model produced nucleus contours. Since the tumor cell detection had superior performance, the detected tumor cells were used to exclude the TIL and tumor nucleus false positive. All computational analyses were executed on a computational server with two CPUs of 22 2.10 GHz cores each, 192 GB memory, and six Nvidia GeForce RTX 2080 Ti GPUs with 11 GB memory each.

We partitioned each WSI into image tiles with a size of 1,024 × 1,024 by pixels to identify tissue regions of high tumor cell density with the DL-based processing pipeline. The top ten image tiles with the highest tumor cell density in each WSI were selected for feature extraction. To generate interpretable models, we chose to select image features of interpretability instead of hidden or intermediate features by machine learning algorithms. Since tumor cells and TILs were reported high correlation with the prognosis or recurrence (37, 38), we extracted three tile-wise features from each image tile, including (1) the tumor cell number, (2) the TIL number, and (3) the tumor cell percentage. Additionally, nuclear grade and TIL number variance were extracted from the ten image tiles collectively. The nuclear grade of each tumor cell was determined by comparing the tumor nuclei size with the adjacent TIL nuclei size. The TIL nuclei size was 304.7 in pixels averaged from representative TILs selected by pathologists. Nuclear grade 1 was defined when the ratio of tumor nucleus size to TIL nucleus size was 1–2.5. Nuclear grade 2 was made when such a ratio was 2.5–3.5. Nuclear grade 3 was made when such a ratio was > 3.5 (Supplementary Figure 2). Tumor cell nuclear grades from the ten image tiles were collected and aggregated to a final nuclear grade by the following rules: (1) if ≥ 10% of the tumor cells had nuclear grade 3, the aggregated nuclear grade was 3; (2) if ≥ 10% of the tumor cells had nuclear grade 2 and rule (1) did not hold, the aggregated nuclear grade was 2; (3) if ≥ 10% of the tumor cells had nuclear grade 1 and neither rule (1) nor (2) held, the final nuclear grade was 1. The image feature of TIL number variance was also computed from the top ten image tiles by cell density as follows:

where V is the TIL number variance; n_i represents the TIL number in the i-th image tile; n¯ is the average TIL number from the ten image tiles. In total, there were 32 image features extracted from each WSI.

A linear regression model was used to correlate with RS. In the regression model, the dependent variable was the RS, while imaging features and Magee features were independent variables. To retain features with high predictive value, we selected features by both domain knowledge and statistical analysis. The independent variables in Magee equations are as follows (39). Magee equation 1 includes Nottingham score, ER and PR H-scores, HER2, tumor size (cm), and Ki67 index; Magee equation 2 includes Nottingham score, ER and PR H-scores, HER2, and tumor size (cm); Magee equation 3 includes ER and PR H-scores, HER2, and Ki67 index. As the feature “HER2” is categorical with two possible values, i.e., “Negative” and “Equivocal,” we used one dummy variable, “HER2_Equivocal,” to represent “HER2” in the regression models. We focused on Magee equation 2 as the Ki-67 index information was missing for more than half samples (195/382, 51.0%) in our datasets. Additionally, the tile-wise features from the first x out of the ten image tiles (x = 1, 2, …, 10), i.e., the tumor cell number, TIL number, and tumor cell percentage, were used jointly. The feature selection was completed in the training set. Various feature combinations were used to construct the linear regression models. The adjusted coefficient of determination R² was used to assess the combinations’ correlation with RS. The feature combination with the highest adjusted R² was selected for the final model.

A total of 7,609 annotated tumor nuclei from 120 testing images and 4,000 annotated TILs from 75 testing images were collected to validate the MRCNN model for tumor nuclei and TIL detection. The trained models correctly detected 6,101 (80.2%) tumor nuclei and 3,304 (82.6%) TILs. Multiple metrics were used for performance assessments, including precision, recall, F1-score, true positive number, false-positive number, and false-negative number. The metrics of precision, recall, and F1-score were defined as follows.

where TP, FP, and FN represent the number of true positive, false-positive, and false-negative samples, respectively. The true positive samples were correctly detected samples. The false-positive samples were cells erroneously detected. Finally, the false-negative samples were missed ground truths from pathologists. The MRCNN models for the tumor nuclei and TIL detection achieved 0.7765 and 0.7171 for the F1-score, 0.7528 and 0.6337 for the precision, and 0.8018 and 0.826 for the recall, respectively.

The Hausdorff distance (HD) was used to measure the tumor nucleus contour concordance between the ground truths from pathologists and predictions using the DL process (Supplementary Figure 3). The metric of intersection over union (IOU) was used to match the ground truth to predicted contours. When IOU was greater than or equal to a cutoff value K, the ground truth and predicted nucleus contours were considered as a matched pair. When there was more than one prediction matching the same ground truth, the prediction with the largest IOU was retained for the match. When one prediction was matched to more than one ground truth, the prediction was assigned to the first matched ground truth. The cutoff value K was set as 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, and 0.9, respectively. We computed the mean equivalent nuclei diameter for each nuclear grade. The mean equivalent diameters for nuclear grades 1, 2, and 3 were 26.30, 34.29, and 48.67 pixels, respectively. We present the mean HD between matched pairs and the ratio of mean HD to the mean equivalent diameter of tumor nuclei in Table 2. Representative cell detection and segmentation results from the DL models are shown in Figure 3.

We detected an overwhelmingly large number of cells in each WSI (Supplementary Table 1). With detection results from image tiles, tumor cell and TIL density distributions were estimated and represented as density maps (Figure 4). The top ten image tiles of each WSI were selected based on the tumor cell density.

Eight variables from the training set included Nottingham grade, ER and PR H-score HER2 status, tumor size (cm), tumor cell number in the densest tile, TIL number variance, and tumor nuclear grade (Table 3). The first five variables were from Magee equation 2, while the last three variables were DL-based image features derived from WSIs. We established a regression model with these selected features from the training set and applied the model to validation sets 1 and 2 for RS correlation.

We divided cases into low, intermediate, and high RS categories with the stratification rules from the TAILORx study (30). The concordances between the RS and our model were 56.10% and 68.0% for validation sets 1 and 2, respectively (Table 4). Additionally, the one-step discordance rates for validation sets 1 and 2 were 39.02% and 48.0%, respectively. The Pearson’s correlation coefficients between the RS and our model were 0.7058 (p-value = 1.32 × 10^–13) and 0.5041 (p-value = 1.15 × 10^–12) for validation sets 1 and 2, respectively. The tumor and TIL density maps from validation sets 1 and 2 are illustrated in Supplementary Figures 4, 5.

The performance of the model correlation with RS was further evaluated by R² and adjusted R² (Table 5). When the image features were integrated with features in Magee equation 2, the adjusted R² value increased from 0.3442 (p-value = 5.17 × 10^–10) to 0.4431 (p-value = 1.32 × 10^–13) in validation set 1 and from 0.2167 (p-value = 6.52 × 10^–12) to 0.2182 (p-value = 1.15 × 10^–12) in validation set 2. Similarly, the R² increased from 0.3846 to 0.4981 in validation set 1 and from 0.2392 to 0.2541 in validation set 2. Additionally, we demonstrated the adjusted R² and R² of the linear regression model that was constructed only with the image features. The resulting adjusted R² and R² are 0.3048 (p-value = 1.61 × 10^–8) and 0.3306 (p-value = 1.61 × 10^–8) for validation set 1 and 0.0139 (p-value = 0.0199) and 0.0309 (p-value = 0.0199) for validation set 2, respectively. It is noted that the image features perform much worse than Magee features in validation set 2. Such performance degradation can be related to the fact that images in validation set 2 were originally scanned at 20 × and later computationally scaled to 40 × magnification. The inconsistency in the original image magnification can contribute to a significant error in the following analyses, leading to a worse prediction result.

To investigate the correlations between Magee and image-derived features, we computed their pair-wise absolute Pearson correlation coefficients. As shown in Figure 5, the largest correlation coefficient of 0.35 was found by the Nottingham score and tumor nuclear grade. Five Magee and image feature pairs present correlation coefficients close to 0.1. All remaining 9 pairs present correlation coefficients less than 0.1. Such weak correlations indicate the complementary prediction value by the image features for RS prediction enhancement.

For further correlation analyses between Magee and image features, we applied the least absolute shrinkage and selection operator (LASSO) regression method to our data and compared the resulting feature coefficients with those in the model trained by Ordinary Least Squares (OLS). The comparison results are presented in Figure 6. As LASSO includes an L1-norm regularizer, it penalizes the excessive feature inclusion and reduces uninformative feature coefficients to zero. From Figure 6, the non-zero feature coefficients from the two models trained by LASSO and OLS present similar values. Coefficients of only three features (i.e., tumor size, HER2, and tumor nuclear grade) were reduced to zero by LASSO. The only removed image feature by LASSO is tumor nuclear grade that presents an absolute Pearson correlation coefficient of 0.35 with the Nottingham score.

We analyzed the cases with discordant risk categories by RS and our model (Table 6). There were totally 54 discordant cases in validation sets 1 and 2. Among these 54 cases, 40 were recommended to have chemotherapy by RS but not by our DL-based model; of these 40 cases, 28 received chemotherapy.

In total, 14 cases were not recommended to have chemotherapy by RS, while our DL-based model did; of these 14 cases, 2 received chemotherapy. The chemotherapy recommendation based on RS and our DL model was determined by the suggested rules from the TAILORx study. Overall, none of these 54 discordant cases developed recurrence regardless of whether received chemotherapy, indicating that the role of chemotherapy in these discordant cases was not clear.