Multiple instance learning-based prediction of programmed death-ligand 1 (PD-L1) expression from hematoxylin and eosin (H&E)-stained histopathological images in breast cancer

Data acquisition and processing

This study was conducted on breast cancer histological samples from three independent cohorts: NF, HY, and XZ (Table 1). The NF cohort consists of 131 patients collected between 2020 and 2024, the HY cohort includes 56 patients collected between 2022 and 2023, and the XZ cohort represents an economically underdeveloped region, with 100 patients collected between 2017 and 2022. The study was approved by the Ethics Committee of People’s Hospital of Xizang Autonomous Region (No. ME-TBHP-22-KJ-041) on July 22, 2023; a waiver of informed consent was granted by the same committee on July 4, 2022. Pathology experts reviewed all available H&E and PD-L1 IHC images and annotated the PD-L1 positive or negative status for all data in these three cohorts. During the annotation process, some patients were excluded from analysis for one of the following reasons: missing H&E or IHC images, non-breast tissue, non-primary diagnosis, or out-of-focus images. The PD-L1 status for the remaining patients was annotated as either negative or positive.

The IHC staining used FDA-approved PD-L1 antibodies (Dako 22C3; Ventana SP263) as companion diagnostics, with PD-L1 (22C3) CPS ≥ 1 and PD-L1 (SP263) TPS ≥ 1 as criteria for PD-L1 positivity. Each slide’s PD-L1 status (positive or negative) was annotated based on the interpretation by clinical pathologists. All slides from the three cohorts were scanned as WSI at 20x magnification using the Teksqray SQS-40P scanner, with a scanner resolution of 0.190 µm/pixel.

In the initial stage of data processing, we implemented a series of rigorous steps to ensure that our model could be trained and validated under optimal conditions.

First, we performed stratified sampling on all case data, dividing it into training (60%), validation (20%), and test cohorts (20%) to ensure that each set represented the characteristics of the overall dataset. Next, to remove interference from background and blurry areas, we applied red, green, blue (RGB) thresholding and canny edge detection to segment the tissue regions in the H&E-stained WSI images, accurately distinguishing white background and blurry areas in the images.

To mitigate the bias from data originating from different centers, we performed normalization as follows. We sliced the WSI images into 512 × 512-pixel patches, with each pixel at a resolution of 0.5 µm. After that, Z-Score normalization was applied to each patch, and the normalization formula is as follows: (1)${\displaystyle y=\frac{x-mean\left(x\right)}{std\left(x\right)}}$

where x is the value of a pixel in a certain channel of the input image, and the mean(⋅) and std(⋅) functions are used to calculated the mean and standard deviation values of the corresponding channels.

The data from different centers may differ in terms of the proportion of positive and negative samples, and the class weight in the binary cross-entropy loss function is used to balance the proportion of positive and negative samples. When calculating the loss, for each sample, the loss is multiplied by the corresponding class weight based on its class. In this way, during back propagation, the model adjusts the parameter update according to the class weight, ensuring that the learning of samples from different categories is balanced.

For biopsy samples, due to the limited sample area, we used a sliding window enhancement algorithm to increase the number of patches. We then applied data augmentations such as rotation and flipping to each patch to effectively expand the variety of patch images.

Finally, to address the imbalance between samples with many patches and those with fewer, we applied random undersampling to the former and random data augmentation to the latter, including operations such as rotation, flipping, scaling, cropping, blurring, Gaussian noise, and salt-and-pepper noise to equalize the patch count for each label. Considering the quadratic complexity of the self-attention mechanism and the large number of patches per WSI image, and given GPU memory limitations (24 GB), we capped the number of patches per label at 12,000.

This series of steps was designed to ensure that our model could fully learn the features within the dataset, thereby achieving strong generalization and improving its performance and accuracy in practical applications.

Local and cross-center validations were performed on the test set, an independent external test set, and the TCGA-TNBC dataset. The models’ prediction performance was compared to confirm their generalizability. Due to the limited number of TNBC patients in the three self-collected datasets, we did not train or validate separately on this group, but instead conducted a statistical analysis of the predicted results. The PD-L1 labels for the TCGA-TNBC dataset were determined based on the study by Jin et al. (2024), where RNA-FPKM values were divided at the 75th percentile; values above this threshold were labeled as PD-L1 positive, and values below this threshold were labeled as negative.

Feature extration

We used the CTransPath model as our feature extraction network to capture a 768-dimensional feature vector for each patch. This model is based on the Swin Transformer architecture, which combines the hierarchical structure of CNNs with the global self-attention module of transformers by computing self-attention within a sliding window. The model starts with three convolutional layers to facilitate local feature extraction and enhance training stability, followed by four Swin Transformer stages. CTransPath was pre-trained on untokenized histopathological images using unsupervised contrastive loss with data from The Cancer Genoma Atlas (TCGA) and the Pathology Artificial Intelligence Platform (PAIP) across multiple organs. This enables the model to act as a collaborative local–global feature extractor, learning generalized feature representations better suited to histopathological image tasks. Each patch’s embedding is stored for the subsequent training process.

Overview of the MIL model

The MIL network architecture for predicting PD-L1 expression in breast cancer histopathological images consists of three main modules: data preprocessing, feature extraction, and the MIL PD-L1 detection network.

First, histopathological images undergo normalized preprocessing to unify image quality and minimize bias due to different scanning conditions. Next, these patches are fed into the CTransPath-based feature extraction network (Wang et al., 2022). The CTransPath network automatically learns deep features within the images, transforming them into high-dimensional vector representations. This high-dimensional vector captures spatial, textural, and structural features in the pathology images, encoding microscopic variations in tumor regions and their surrounding environment.

The high-dimensional features are then aggregated and analyzed by a weakly supervised PD-L1 detection network based on MIL. The MIL framework learns the relationship between multiple patch features and PD-L1 expression, effectively handling information inconsistencies across patches, thereby enabling accurate PD-L1 expression prediction under weak supervision for each histopathological image. Finally, the model produces a PD-L1 prediction based on the aggregated features of multiple patches and outputs a confidence score to assess the reliability of the result.

Heatmap visualization of interpretable features

The heatmap is generated by extracting the image patch features in the WSI and combining the model-generated attention scores to map to the original coordinates of the WSI. The scores can be normalized to percentiles for comparison, resulting in a visual heat map that covers the entire WSI and visually shows the model focus areas and their importance. The MIL-based approach enables the generation of highly interpretable heatmaps. During model training, pixel-level annotations are not required, yet the model can intuitively display the relative contribution and importance of each tissue area to its predictions, facilitating clinical analysis. We used a feature aggregation model that assigns an attention score to each patch.

In the heatmap, higher attention values (shown in red) indicate areas that contribute more significantly to the model’s prediction, while lower attention values (shown in blue) represent less impactful regions. By mapping these attention values across the whole image, the resulting heatmap clearly highlights the regions the model focuses on most.

This heatmap is a valuable tool for helping clinicians understand the model’s decision-making process. Red areas suggest regions the model considers more likely to be associated with lesions or critical pathological features, whereas blue areas are seen as less relevant to the prediction. This visualization not only enhances model interpretability but also provides meaningful guidance for clinical diagnosis.

Statistical analysis

The model classification results are visualized as a confusion matrix. According to the confusion matrix, we use area under curve (AUC), accuracy (ACC), precision (PRE), specificity (SPE), sensitivity (SEN) and F1-score as our statistical metrics to quantify the performance of the method. AUC with 95% confidence intervals (CI) were calculated to evaluate the stability and accuracy of these models. The specific formulas are shown as Eqs. (2)–(6). (2)${\displaystyle ACC=\frac{TP+TN}{TP+TN+FP+FN}}$ (3)${\displaystyle SEN=\frac{TP}{TP+FN}}$ (4)${\displaystyle SPE=\frac{TN}{TN+FP}}$ (5)${\displaystyle PRE=\frac{TP}{TP+FP}}$ (6)${\displaystyle F\text{-}1=\frac{2\times SEN\times PRE}{SEN+PRE}}$

where TP = true positive, TN = true negative, FP = false positive, and FN = false negative.

Experimental setup

Experiments were conducted using Python 3.8 and PyTorch on a server equipped with a single NVIDIA 4090 GPU and an Intel Core i9-14900K CPU. This hardware configuration and software environment ensured efficiency and consistency during the model training process and provided high-performance computing resources to support the complex computational requirements of the models.

The Adam optimizer was used for all three models training. The weight decay parameter of each model was 1e–2 for AttentionMIL, 1e–4 for CLAM and 1e–5 for TransMIL. While the learning rate was 5e–4 for AttentionMIL, 2e–4 for CLAM and 5e–4 for TransMIL.

All models were trained for a total of 30 epochs. Due to GPU memory limitations, the batch size was set to 8. In each epoch, we employed the MultiStepLR learning rate adjustment strategy, reducing the learning rate to 0.1 times the original value at the 15th and 25th epochs to enhance the convergence and stability of the models. This optimized setup was then applied to evaluate the models on the test set. In order to guarantee the accuracy of the experimental results, all experimental results are the average of three runs.

Experimental results

We developed three deep learning models: AttentionMIL, CLAM, and TransMIL. We tested the ability of these three models to predict PD-L1 expression directly from pathology images using the NF and HY breast cancer datasets. Due to the limited size of the datasets, we divided them into 60% training, 20% validation, and 20% testing sets. The results showed that the TransMIL model achieved the highest AUC of 0.833 on the NF and HY test sets, outperforming CLAM (0.815) and AttentionMIL (0.806) (Figs. 2A–2D). The scatter plots visualized the classification of the results of each model on the test sets (Figs. 3A–3C). To further evaluate the stability and accuracy of AUC, the 95% CI of the three models in the two test sets were calculated (Table 2). The results showed that the 95% CI of the three models in the two test sets were all around 0.1, and the classification performance of the models was relatively stable. However, the AUC 95% CI of the three models had a large overlap, and the differences between the models were not significant. Therefore, we further calculated the accuracy, precision, sensitivity, specificity and F1-score of each model (Table 3). In the NF and HY test sets, AttentionMIL and CLAM performed similarly, but both had slightly lower accuracy and precision compared to TransMIL, while their sensitivity lagged behind TransMIL, indicating that TransMIL performed best in detecting positive samples (Fig. 3G).

Given the significant differences in images between surgical and biopsy samples, we used the best-performing TransMIL model to separately analyze the prediction results for both surgical and biopsy samples. The results showed that the mean AUC for surgical samples was 0.856, while the AUC for biopsy samples was 0.762 (Fig. 4B). This indicates that the model performs better, or at least as well, for surgical samples. Together, these results demonstrate that our models can predict PD-L1 expression levels from H&E images, with better performance for surgical samples than for biopsy samples.

Generalization performance

To further validate the generalization ability of the models, we used an independent external validation set, the XZ breast cancer biopsy data, for testing. The data from the XZ cohort exhibits distinct characteristics compared to those of NF and HY datasets. The TransMIL model achieved an AUC of 0.799 on the XZ test set, showing a decrease of approximately 4.1% compared to the NF and HY test sets (Fig. 4A). CLAM and AttentionMIL models also showed a similar decline in performance on the XZ test set, with AUC values of 0.786 and 0.778, respectively (Figs. 2E–2H). The scatter plots were visually displayed the classification of the results of each model on the XZ test sets (Figs. 3D–3F). We also calculated the accuracy, precision, sensitivity, and specificity for each model (Fig. 3G). In terms of accuracy, precision, and specificity, TransMIL outperformed both AttentionMIL and CLAM, demonstrating superior classification ability. However, AttentionMIL achieved the highest sensitivity, indicating its superior performance in detecting positive samples in the XZ test set.

When the best-performing TransMIL model was applied to analyze the predictions for surgical and biopsy samples in the XZ test set, the results showed an AUC of 0.822 for surgical samples, compared to 0.750 for biopsy samples (Fig. 4C). These results indicate that all three models demonstrate good generalization ability on the independent external test set XZ, with TransMIL still performing the best, and surgical samples outperforming biopsy samples.

Studies have shown that using extensive RNA sequencing data can predict PD-L1 expression in H&E-stained slides. This suggests that RNA sequencing results have a threshold that allows H&E pathological images to be distinguished. Given the limited number of TNBC cases in this study (only 19 cases), the TransMIL model achieved an AUC of 0.750, while AttentionMIL and CLAM performed lower, with AUC values of 0.716 and 0.682, respectively (Figs. 5A–5D). Therefore, we aimed to use TCGA TNBC data to validate the consistency between our models and RNA-based prediction models.

Based on previous study (Jin et al., 2024), we selected TNBC cases from the TCGA cohort and set the RNA FPKM value threshold as the 75th percentile, which was 1.9355. Values above this threshold were classified as positive, while those below were classified as negative. In the TCGA TNBC dataset, the TransMIL model achieved an AUC of 0.721, while CLAM and AttentionMIL performed lower, with AUCs of 0.639 and 0.633, respectively (Figs. 5E–5H). The prediction performance of all three models in the TCGA-TNBC dataset was lower compared to the NF and XZ datasets (Figs. 5A, 5E).

A scatter plot of the prediction results from the best-performing TransMIL model revealed that the model had very few false negatives but a relatively high number of false positives (Fig. 6A). We also calculated accuracy, precision, sensitivity, and specificity for the predictions of the NF and XZ_TNBC and TCGA-TNBC datasets. The results showed that all three models had high sensitivity, even reaching 1, but their specificity was relatively low, especially for AttentionMIL and CLAM in the TCGA-TNBC dataset (Fig. 6B).

Although the high sensitivity (close to 1) indicates that the models performed well in detecting positive cases, the low specificity, particularly for TransMIL which performed below 0.5, reflects a high false-positive rate. While this may reduce the chance of missing positive cases in clinical practice, it could also lead to an increased economic burden and time costs for a significant number of false-positive patients. Overall, we believe that the TransMIL model still performed well in predicting PD-L1 expression in TCGA-TNBC H&E slides, but further optimization and adjustments are needed for broader clinical application.

Comparison of performance differences between models

Multiple studies have proven that PD-L1 is an important biomarker for immune therapy response, with significant implications across various types of cancer (Chen & Mellman, 2013; Upadhaya et al., 2022; Yi et al., 2022). The FDA-approved PD-1/PD-L1 immunotherapies currently include treatments for lung cancer, colorectal cancer, TNBC, and other malignancies (Lee et al., 2020; Nobin et al., 2024; Savic et al., 2019). In the rapidly developing field of digital pathology, many tools for predicting molecular markers from tumor pathology slides have been developed (Baxi et al., 2022; Bera et al., 2019; Shafi & Parwani, 2023). In this study, we trained and made prediction using three models: AttentionMIL, CLAM, and TransMIL.

Whether on the NF and HY test sets or the independent XZ test set, the ACC and F1-scores of TransMIL were higher than those of CLAM and AttentionMIL, indicating its overall better classification performance. Sensitivity and specificity measure how well positive and negative examples were identified, respectively. From the experimental results, in the NF and HY test set, TransMIL had the highest sensitivity, but in, but in the XZ test set, TransMIL had the same sensitivity as CLAM and slightly lower sensitivity than AttentionMIL. In terms of specificity, the three models performed consistently on the NF and HY test set, but TransMIL was significantly higher than the other two models on the XZ test set. The sensitivity was high, indicating that the model performs better in detecting positive samples, but AttentionMIL and CLAM had significantly lower specificity, reflecting that their good detection of positive samples might be at the expense of detecting negative samples, and the model performance had a significant deviation. The accuracy also clearly reflects this feature, and in the NF and HY and XZ test sets, especially in the latter, TransMIL showed better results, demonstrating its better generalization ability.

For such differences, researchers believe that, on the one hand, it may be related to the characteristics of the Transformer architecture. The Transformer architecture in TransMIL better utilizes the feature relationships between different patches, capturing global information about biomarkers within tissue samples, making it more suitable for the label generation methods used in NF and HY. On the other hand, the Transformer-based long-distance feature aggregation mechanism helps the model better adapt to the data characteristics of Center 3 to a certain extent. In contrast, CLAMMIL’s instance selection mechanism achieves a slightly lower AUC on Center 3, indicating that the instance selection strategy may be sensitive to more heterogeneous data when sample labels are generated in a way that deviates from the training data. On the other hand, AttentionMIL only uses attention weights, which makes it difficult to fully adapt to changes in data distribution and label generation methods, so it has relatively weak performance in cross-center testing.

Meanwhile, compared to the pre-trained NF and HY datasets, all models showed decreased performance on the independent external XZ test set. This may be due to differences in how labels were obtained. In NF and HY, PD-L1 labels were acquired individually from tissue samples using conventional full-slide immunohistochemistry methods, whereas XZ used a tissue microarray (TMA) method for immunohistochemical detection. The tissue microarray method requires high selectivity and tissue homogeneity, which may not fully represent the entire tissue sample, potentially leading to information loss. This could explain the performance decline of the models on the XZ test set. In addition, despite data normalization, there are differences in staining instruments and reagents across different centers, differences in data storage formats, and sample bias caused by variations in race and living environment among datasets, which may lead to the performance degradation of each model on the XZ test set.

Comparison of differences between datasets

This study involves data from three centers, primarily using formalin-fixed paraffin-embedded (FFPE) samples, including both biopsy and surgical samples. This provides a solid sample foundation for exploring the stability of model performance across multiple centers and sample types. As an independent external validation set, the XZ dataset comes from the Xizang Autonomous Region in China, which is relatively underdeveloped compared to the economically advanced coastal regions. Medical resources in the region are limited, and maintaining consistent quality in pathological slides is challenging. Additionally, the cohort is predominantly composed of the local Tibetan ethnic group. Given the unique high-altitude environment, cold climate, and distinct lifestyle, dietary habits, and other factors, the XZ cohort holds significant importance.

The experimental results show that the overall predictive AUC reaches around 0.800, demonstrating good generalization, with surgical samples showing better predictive accuracy than biopsy samples. On the one hand, surgical samples contain more information, including not only tumor regions but also tumor invasion edges, tumor-associated lymphocyte infiltration, and peritumoral stroma. On the other hand, biopsy samples are more technically demanding for the operator, requiring accuracy in the biopsy site, and are limited in size, potentially missing important surrounding tissue information. These factors may contribute to the observed performance differences between surgical and biopsy samples. Analyzing the locations of the three centers used in this study reveals that Center NF and Center HY are both located in the economically developed southern coastal regions, where sample collection is relatively more abundant, and the quality control during processes like sampling, preparation, and staining is stricter, resulting in more consistent samples. In contrast, Xizang is located in the inland, economically underdeveloped region with limited medical resources, and quality control in pathological slide preparation is relatively weaker. Despite image preprocessing steps like staining normalization, these disparities may still lead to differences in model performance.

Visual heatmaps showing interpretability features

By visualizing heatmaps, we further explored the connection between model predictions and histological features. In molecular marker prediction tasks, heatmaps can highlight tissue regions that are highly correlated with specific molecular markers (e.g., PD-L1 in this study) and analyze the relationship between tissue morphology and molecular marker expression, enabling pathologists to visualize the basis of the model’s predictions and helping to understand how pathological features can influence changes at the molecular level. The attention heat map reveals the tissue structures that the model focuses on when making predictions, enhancing the transparency and interpretability of the model. In addition to regions of interest, regions of uncertainty and their confidence levels are also visualized, thus assisting pathologists in their review and reducing the risk of misclassification.

Previous studies have shown that features such as high nuclear heterogeneity, tumor necrosis, lymphovascular invasion, inflammatory cell infiltration, stromal reaction, tumor necrosis, and high-level tissue patterns are important in the diagnostic nuclear prognostic assessment of solid tumors and often correlate with poor tumor prognosis (Arpinati & Scherz-Shouval, 2023; Bilotta, Antignani & Fitzgerald, 2022; Stanton & Disis, 2016). In this study, the best predictive TransMIL model was subjected to heatmap visualization to observe the histological features of the focal patches to obtain clinical interpretability (Figs. 7A–7D). In the heatmap visualization results of TransMIL, the patches with high predictive scores exhibited features such as tumor cell nuclei with large and pronounced nucleoli, lymphocyte infiltration, and high-grade histological patterns (solid or micropapillary structures), which were consistent with previous study results. The features of large nuclei and pronounced nucleoli of tumor cells tend to be associated with features such as high proliferative activity, high malignancy, and genetic mutations, which in breast cancer often also suggests a possible association with chemotherapy tolerance. Tumor infiltrating lymphocytes play an important role in the prognostic assessment of breast cancer, and their effect on tumor is associated with distribution location, cell subtype, and tumor subtype. The model focuses on the characteristics of lymphocyte infiltration, which, while suggesting the site of abnormality, can also help clinicians to predict the response to treatment in conjunction with tumor typing. In addition, the high-level histological pattern focused on the heatmap was classified as the relatively worst group in terms of prognosis among the high, medium and low histological patterns classified by the World Health Organization (WHO) as recently as 2015 (Moreira et al., 2020; Travis et al., 2015). All of the above histological features suggest that the TransMIL model, which was the most predictive in this study, facilitates focusing on the lesion site and also allows for the acquisition of histological features associated with a poor prognosis, such as high-level histological patterns. The acquisition of such histological features allows the model to have more accurate predictive performance while providing clinicians with interpretable features for diagnosis and prognosis.

Model performance in TNBC subgroup

Currently, clinically approved PD-L1 testing and corresponding treatments are primarily targeted at TNBC patients in breast cancer, while first-line treatments for non-TNBC patients still involve hormone therapy or Her-2 targeted therapy (Emens & Loi, 2023; Prat et al., 2015). To explore the predictive effectiveness of our model for TNBC patients, we classified patients based on ER, PR, and HER-2 information, and made predictions for a total of 19 patients from the NF and XZ centers. The results showed that the TransMIL model performed the best. The confusion matrix clearly indicates that all three models tend to predict TNBC patients as positive. In clinical applications, this would significantly reduce the rate of missed diagnoses. Since our TNBC sample size is very limited, it is difficult to determine how the model would perform on a larger scale. However, studies have shown that RNA sequencing can predict PD-L1 expression in H&E slides. Therefore, we also made predictions on the publicly available TCGA-TNBC dataset, and the results showed an AUC of 0.721. These results suggest that, not only in the TNBC dataset but also across all datasets, our model tends to identify patients as positive. We hypothesize the following reasons for this. First, the distribution of the test set differs from that of the training set, where the training set has more positive samples than negative ones, while both the independent external test set (XZ) and the TCGA-TNBC dataset have more negative samples than positive ones. The uneven sample distribution between the training and test sets, along with these discrepancies, may cause the model to lean towards predicting positivity. Secondly, although data from three centers were included, the sample size from each center is very limited, which may restrict the model’s performance. Additionally, models trained on RNA sequencing results have shown an AUC of only 0.540 when predicting PD-L1 expression in FFPE pathological slides from TNBC. Even with better performance using fresh frozen tissue slides, the AUC is only 0.640. This suggests that while there is a morphological-transcriptomic link in PD-L1 expression, its separability in TNBC is generally limited. In the future, it is hoped that using more integrated morphological-transcriptomic data for training will improve the model’s predictive performance.

Innovations and limitations

In summary, this study combines multi-instance learning with the XZ dataset for integrated innovation. The MIL approach allows the AI to automatically learn key regions from WSI, freeing it from the dependence on pixel-level annotation, and improves the model’s generalization ability while effectively reducing the workload of manual annotation. The XZ dataset is from the Tibetan region, where medical resources are comparatively limited and the cohort is predominantly composed of ethnic minorities. This results in disparities among cohort populations and in the level of preparation and reading of pathology slides. The model’s reliable cross-center capability enables the implementation of training using data that is readily available in developed regions and subsequent prediction of data from less developed regions. This is important for achieving the efficient allocation of healthcare resources in regions with unbalanced economic development.

However, this paper also has some limitations. The first is the issue of data imbalance. Due to the limited size of the dataset, PD-L1 expression is unevenly distributed across the datasets, especially the TNBC subgroup, which shows opposite distribution characteristics in the training set and the independent external test set. Secondly, three local datasets were used in this study, with data originating from different pathology laboratories, each with differences in tissue processing methods, staining techniques, and technical practices of personnel. Despite data normalization, heterogeneity still existed among the data. Moreover, it is clear from the results that the performance of all three models on the XZ dataset is degraded, suggesting that race-specific and region-specific datasets may have an impact on the generalizability of the models. Finally, although the model achieved AUC up to about 0.800 in the local dataset, it still lacks clinical validation tests in a really patient population. In the future, we will increase the data size, address data imbalance and improve the normalization process, increase the model generalizability and deploy it in the clinic, applying and improving the model in real environments with a view to its early use.