A multi-class classification algorithm based on hematoxylin-eosin staining for neoadjuvant therapy in rectal cancer: a retrospective study

Patient selection

A total of 94 patients who received NAT in Chongqing University Cancer Hospital (CUCH) from 2011 to 2022 were retrospectively enrolled. The criteria for inclusion comprised: (1) Diagnosis of primary rectal cancer through pathological biopsy before NAT; (2) receiving of the NAT regimen based on Capecitabine or 5-FU combined with Oxaliplatin; (3) availability of preoperative pathological biopsy sections; (4) have underwent surgery after NAT and have TRG evaluation results. The criteria for exclusion comprised: (1) Other treatments before NAT; (2) lack of preoperative pathological biopsy section; (3) poor quality of H&E-stained images. This study was approved by the Ethics Committee of Chongqing University Cancer Hospital (Ethics number: CZLS2022185-A).

Pathological assessment of response and data collection

Specimens were prepared onto glass slides and stained with hematoxylin & eosin (H&E). TRG grades were evaluated by two experienced pathologists, and samples with discordant scores were reevaluated until an agreement was reached. In addition, clinical data including the age at diagnosis, gender, tumor/node(T/N) stage, CEA and CA199 levels before treatment, and treatment regime were collected. Due to the institutional policy, the pathological image data are not publicly available.

Image preprocessing

H&E-stained slides were produced at 40× magnification by the KFBio KF-PRO-005-HI digital scanner. There were 133 slides collected. Whole slide images (WSIs) were annotated using QuPath open-source software (version 0.3.2). The annotation of tumor tissues in the WSIs was performed by two pathologists. One pathologist was responsible for delineating, and the other for examining. To reduce the computing consumption, the WSIs and their corresponding annotations were split into overlapping 512 × 512 pixels tiles at 0.238 µm/pixel resolution (equivalent to 40× magnification) and 1 µm/pixel resolution (equivalent to 10× magnification) with overlapping segmentation, which stride was0.5. Tiles with tissues less than 50% were removed manually.

All tiles were color normalized with the structure-preserving color normalization (SPCN) method (Vahadane et al., 2016). Different types of data augmentation were applied, including random vertical flip (0.5), random horizontal flip (0.5), and brightness/contrast/saturation/hue changes (1).

Construction and visualization of response classifier

The WSIs from CUCH were split at a ratio of approximately 4:1, with 113 WSIs in the training set and 20 WSIs in the test set. The training set consisted of 166,919 tumor ROI tiles at 40× magnification. There were 14 slides (27,598 tiles) from TRG0, 70 slides from TRG1/2 (98,102 tiles), and 28 slides (41,219 tiles) from TRG3. We included multiple slides per patient, a patient had one or two slides. After a series of data augmentation techniques, the augmented training set increased to three times the original training set. The test set was seven slides (15,506 tiles) from TRG0, seven slide (13,407 tiles) from TRG1/2, and six slides (10,263 tiles) from TRG3. For feature extracting, we utilized a deep residual network (Resnet101) without the final softmax classification layer and average pooling layer, pretrained with ImageNet. After the 2048-dimension feature vector were extracted, multiple full connected layers were used to non-linear mapping. At the end, a classification layer was added. In the training process, the tiles with 512 × 512 pixels were re-sized to 256 × 256 pixels and then input into the network. Among them, the patient’s treatment responses were considered as the target label. We optimized the response classifier using a five-fold cross-validation approach, which involves randomly dividing the training set into five balanced subsets. Based on TRG classes, each fold was divided into 4:1 from different TRG class to extract training samples and validation samples. The hyperparameters include a batch size of 32, an initial learning rate of 0.001, and total epochs of 20 (decreasing learning rate by 0.1 every five epochs). In the post-processing stage, the average score of the tumor tiles of each WSI was taken as the prediction probability of the WSIs, which was regarded as the criterion to judge the efficacy of NAT. The response classifier structure was implemented in Python, using Pytorch backends (Python 3.9.7, Pytorch 1.11.0, SKLeran 0.24.2). We used the same network structure, batch size, and learning rate to train and verify the model at 10× magnification. The initialized parameters with pre-training weights in ImageNet.

To analyze features used by the classifier to make its decision, we performed Smooth Grad-CAM (class activation mapping) (Omeiza et al., 2019) to identify regions within each tile that the neural network uses to generate predictions. Then, we used cell detection and classification functions in Qupath to identify three types of cells: tumor cells, tumor-infiltrating lymphocytes (TILs), and stroma (Acs et al., 2019). TILs% calculated as follows: ${\displaystyle TILs\text{\%}=\frac{TILs}{TILs+Tumorcells}\times 100\text{\%}.}$

Statistical analysis

We evaluated models at the tile and slide levels. The neural network yielded a probability value for every class of interest in each tile. Then we averaged the probabilities of each tile from the patients’ slides to assign a final probability to each patient. The probabilities calculated by the neural network were used to generate the area under the curve (AUC) for the receiver operating characteristic (ROC) curve. Using macro-AUC and micro-AUC as metrics, we compared the performance of the models at different magnifications. We also used Precision, Sensitivity, Specificity, and the multi-class Matthews’ correlation coefficients (MCC) to assess the performance of models. The statistical analyses were performed using SKLeran 0.24.2 in Python (scikit-learn.).

Workflow for predicting NAT response with a multi-class classification algorithm

A total of 94 patients were included in this study based on the inclusion/ exclusion criteria. Table 1 summarizes the characteristics of TRGs of the cohort. The workflow of our studies is shown in Fig. 1. Firstly, we collected and preprocessed the HE images from the patients, including scanning, delineating and segmenting the ROIs. In this step, a total of 166, 919 tiles were obtained. Then we normalized the color of tiles by the SPCN (Vahadane et al., 2016). Secondly, we input the tiles into the neural network in the form of tiles at 10× and 40× magnification, and extract the features through ResNet-101 to deduce the probabilities and predict the treatment responses: TRG0, TRG1/2, or TRG3. Finally, we constructed a visual heat map of tiles using CAM to explore the medical interpretability of the model.

Model performance evaluation under multi-class classification tasks on different scale

To test the model performance on different fields of view, HE tiles were preprocessed with two magnifications, including that of 40× and 10×. We first developed the response classifier at 40× magnification. At the tile level, the macro-AUC of the three-class classification model was 0.82, and the micro-AUC was 0.86 on the cross-validation of the training set, and the macro-AUC of 0.72 and the micro-AUC of 0.72 on the test set, as shown in Figs. 2A and 2E. The model revealed a stable performance with the macro-AUC of 0.97 and the micro-AUC of 0.96 at the slide level, as shown in Fig. 2B. In order to compare the performance on different scales, with this dataset we trained and tested the classifier at 10× magnification. The AUCs of the tile and slide levels are shown in Figs. 2C and 2D. The macro-AUC and micro-AUC were 0.75, 0.78 and 0.79, 0.86, respectively. And the result of the test set and the variability on the 5-folds cross validation are shown in Figs. 2E–2H, Tables 2, and 3.

After comparing the overall performance of models, we determined the specific classification performance of the model for each group shown in Fig. 3. Specifically, the performance of the model at the 40× magnification is better than at the 10× magnification, which indicates that the local cellular features extracted from the model were more conducive to the differentiation of curative effects than the global tissue features. The detailed AUCs under tile and slide level are provided in Table 4.

Medical interpretability analysis

To explore the relationship between histopathological features and the efficacy of NAT, we developed a post-hoc explanation of a trained model at the 40× and 10× magnifications, based on the CAM algorithm. The model generated a heat map for each tile, and the areas with large contribution from the model are displayed in red. As shown in Figs. 4A and 4B, overlaying the original tile images with the image of the CAM analyses revealed that tumor tissues including tumor cells and TILs seem to contribute the most in our model, while stroma contributes less to this model. Among them, TILs play an important role in the decision to classify responses, either with 40× or 10× magnifications. Specifically, tiles predicted with TRG0 appeared to be denser in the number of TILs than that with TRG1/2 or TRG3. As shown in Fig. 4E, it was observed that the proportion of TILs was significantly associated with the treatment response. The median TIL % in the TRG0 group was significantly higher than that in the other two groups (the TRG0 group was 35.1%, the TRG1/2 group 24.7%, and TRG3 group 16.9%).

Correlation analysis of TRG prediction

As shown in Table 5, at the 40× magnification, the MCC in the prediction of TRG0,1/2 and 3 at the tile level were 0.499, 0.501 and 0.425, respectively, while the MCC in the prediction of TRG0, 1/2 and 3 at the slide level were 0.617, 0.734, and 0.683, respectively. At the 10×  magnification, the MCC in the prediction of TRG0,1/2 and 3 at the tile level were 0.623, 0.233 and 0.221, respectively, while the MCC in the prediction of TRG0, 1/2 and 3 at the slide level were 0.690, 0.316, and 0.298, respectively. According to the MCC results, TRG1/2 and TRG3 were weakly differentiated at low magnification. We also found that the proportion of TILs in TRG1/2 and TRG3 was lower than that in TRG0, and this microstructure could not be clearly displayed at low magnification if the proportion was low. Therefore, this may be the reason for the low MCC in pathological images with low magnification. Meanwhile, according to the MCC analysis results, we found that TRG0 had better predictability at low magnification, while TRG1/2 and TRG3 hah relatively better predictability at high magnification. Therefore, combination of low magnification and high magnification pathological images to construct a multi-scale joint prediction model might improve the accuracy and robustness of TRG classification prediction for NAT.