Breast cancer homologous recombination deficiency prediction from pathological images with a sufficient and representative Transformer

Abstract

Homologous recombination deficiency (HRD) has been recognized as a key biomarker for poly-ADP ribose polymerase inhibitors (PARPi) and platinum-based chemotherapy in breast cancer (BC). HRD prediction typically relies on molecular biology assays, which have a high turnaround time, and cost. In contrast, tissue sections stained with hematoxylin and eosin (H&E) are ubiquitously available. However, current HRD prediction methods that utilize pathological images are usually based on attention-based multiple instance learning, which is ineffective for modeling the global context of whole slide images (WSIs). To address this challenge, we propose a Sufficient and Representative Transformer (SuRe-Transformer) for WSI-based prediction of HRD. Experimental results demonstrate the superior performance of SuRe-Transformer in predicting HRD status compared to state-of-the-art methods, achieving an AUROC of 0.887 ± 0.034. Furthermore, SuRe-Transformer demonstrates generalizability across multiple external patient cohorts and achieves state-of-the-art performance in predicting several gene mutation biomarkers from BC WSIs.

Introduction

Breast cancer (BC) has the highest incidence rate among women globally and continues to be a leading cause of cancer-related mortality in women¹; advancing the treatment of BC is of paramount importance. Homologous recombination repair (HRR) is a high fidelity repair pathway of double-strand DNA breaks. Homologous recombination deficiency (HRD) has been recognized as a key biomarker for poly-ADP ribose polymerase inhibitors (PARPi) in BC. Several methods have been developed to detect HRD, including testing for mutations in HRR-related genes², analyzing genomic instability³ and mutational signatures⁴, and conducting functional assays⁵. However, standard diagnostic tests for detecting HRD via molecular profiling, pose economic challenges, particularly in low-resource settings⁶. Consequently, there is an urgent need to develop affordable and available methods for HRD prediction.

Hematoxylin and eosin (H&E) stained pathological images are ubiquitously available in clinical practice. Recent work by Lazard et al.⁷ indicates that morphological features, including high tumor cell density, conspicuous nucleoli, tissue necrosis, distinctive laminated fibrosis, and tumor infiltration, serve as the key indicators for differentiating HRD and homologous recombination proficient (HRP). Therefore, we decide to predict HRD status from WSIs stained with H&E. WSIs typically exhibit extremely high resolution and are often stored in a multi-resolution pyramid structure⁸. Consequently, it is generally considered ineffective to directly apply deep learning methods to the entire WSI for training⁹. To address this issue, we select several patches from WSI and aggregate these patches as WSI-based feature representations for HRD status prediction. While several approaches^10,11 randomly select patches from WSIs, and others utilize color-based methods¹² or patch-level classifier predictions¹³, these methods often fail to ensure the representativeness of the selected patches. Therefore, a good patch selection method, especially one that ensures the representativeness of the selected patches, is necessary for accurate HRD status prediction.

In the evolving landscape of bioinformatics, machine learning approaches are progressively demonstrating transformative potential^{14,15,16,17,18,19,20,21}. This trend is particularly evident in computational pathology, where recent advances in deep learning²² have demonstrated accurate, reliable, and reproducible performance across various pathological diagnosis tasks, including the capacity to identify patterns not typically recognized by human experts. Schirris et al.²³ introduce DeepSMILE, which uses a two-stage multiple instance learning (MIL) that combines self-supervised learning (SSL) for convolutional neural network (CNN) pre-training to extract features and variability-aware multiple instance learning (VarMIL) for feature aggregation to predict HRD status. To mitigate a bias among computational pathology workflows, Lazard et al.⁷ introduce a strategic sampling approach. They utilize Momentum Contrast (MoCo) for feature extraction and aggregate these features based on MIL to predict HRD status from tissue slides and identify related morphological phenotypes associated with HRD. Bergstrom et al.²⁴ propose DeepHRD, which is based on ResNet18 and MIL, to predict the HRD status from BC WSIs.

Currently, although progress has been made in MIL-based HRD prediction methods using pathological images, these methods are ineffective in capturing the global context of WSIs. However, a new type of model architecture has emerged to replace the attention-based MIL method: Transformers. Originally developed for natural language processing (NLP), Transformers²⁵ are robust tools for sequence processing, effectively harnessing the potential of large-scale data. Following the success of Transformer in NLP, Transformer has recently achieved remarkable success in various visual recognition tasks, primarily due to its ability to effectively capture long-range dependencies in sequences through attention mechanisms²⁶. Dosovitskiy et al.²⁷ proposed ViT, which directly applies a Transformer architecture to non-overlapping image patches for image classification. However, vanilla Transformers are constrained by computational complexity and can only handle shorter sequences of fewer than 1000 tokens²⁸. To address these challenges, several Transformer-based methods for WSI analysis have adopted linear attention mechanisms^29,30, including the Nyströmformer²⁸ and Performer³¹, which replace the traditional softmax attention to reduce memory usage. Additionally, several studies have adopted sparse attention³² and deformable attention³³ to reduce computational costs and memory requirements. However, the above Transformer-based methods often encounter overfitting issues when dealing with limited pathological images due to the nature of self-attention.

In the present study, we propose a novel Transformer-based framework, the Sufficient and Representative Transformer (SuRe-Transformer), for predicting HRD status from BC WSIs, as illustrated in Fig. 1. Altogether, we demonstrate that SuRe-Transformer addresses several key challenges in applying the Transformer to WSI through the following methods: (1) Instead of randomly selecting patches from WSI, we employ cluster-size-weighted sampling to ensure the representativeness of the selected patches. Additionally, we provide a mathematical proof for the effectiveness of this patch selection method. (2) We propose a radial decay self-attention (RDSA) to extend the input sequence length in Transformer, ensuring that a sufficient number of patches are considered. (3) We adapt an unsupervised feature extractor, pretrained using DINO on the largest BC WSI dataset we curated, to extract patch embeddings, enabling better feature representation for clustering and Transformer. Experiments conducted on BC from The Cancer Genome Atlas (TCGA) dataset demonstrate the superior performance of SuRe-Transformer in predicting HRD status compared to state-of-the-art methods. Furthermore, predictions from SuRe-Transformer generalize effectively to an external BC cohort. Additionally, our proposed framework can capture the general morphological features related to HRD, as demonstrated by an HRD classifier trained on TCGA BC that is able to predict HRD in TCGA ovarian and pancreas cancers. Moreover, our proposed framework is a general WSI representation learning framework that can predict various gene mutation biomarkers, including TP53, GATA3, MAP3K1, ERBB2, BRCA1, and PIK3CA, from BC WSIs.

Fig. 1: Workflow overview with cohort & pretraining, patch selection and mode architecture.

a Overview of the 6 cohorts of BC, which are used in the study for pre-training feature extractors using self-supervised learning. Created with BioRender.com. b The preprocessing workflow for WSIs consists of several steps: WSI digitization, tissue segmentation, and tessellation, patch embedding based on self-supervised learning, and patch selection through clustering and cluster-size-weighted sampling. c The illustration of the SuRe-Transformer and network details of the SuRe-Transformer encoder.

Results

The experiment in this study is divided into four parts. First, we use TCGA BC dataset to compare SuRe-Transformer with state-of-the-art methods in evaluating the HRD status. Second, we conduct three external tests to compare SuRe-Transformer with DeepSMILE, as shown in Fig. 2a. Third, we evaluate the prediction performance of SuRe-Transformer in predicting six gene mutation biomarkers on TCGA BC dataset, as illustrated in Fig. 2b. Finally, we provide a visualization analysis to illustrate the results of patch clustering and conduct an ablation study to assess the contribution of each component of SuRe-Transformer.

Fig. 2: Flowchart of all experiments used for evaluating SuRe-Transformer performance.

a The area under the receiver operating characteristic curve provides an overall measure of SuRe-Transformer accuracy in predicting homologous recombination deficiency, evaluated using held-out test sets comprising different cancer types within a single data center (TCGA-PAAD, TCGA-OV) as well as data from distinct centers for the same cancer type (ChosenMed-BRCA). b The area under the receiver operating characteristic curve provides an overall assessment of SuRe-Transformer accuracy in predicting gene mutation biomarkers-including TP53, PIK3CA, MAP3K1, GATA3, ERBB2, and BRCA1-from breast cancer whole slide images. Created with BioRender.com.

Dataset description

We acquire Formalin-Fixed, Paraffin-Embedded (FFPE) WSIs for BC from TCGA. The dataset, collected from multiple medical centers across the United States, is available for download at https://portal.gdc.cancer.gov/. A description of the patient cohort is shown in detail in supplemental material (Supplementary Tables 1–7). We retrieve genomic DNA Damage Response (DDR) related labels, including HRD Score, from Knijnenburg et al.³⁴. The HRD score is calculated based on several criteria: the number of subchromosomal regions with allelic imbalance extending to the telomere³⁵, the number of chromosomal breakpoints between adjacent regions of at least 10 megabases³⁶, and the number of regions with a loss of heterozygosity event of intermediate size³⁷. The methods for labeling HRD status are inspired by binary (mHRD) and ternary (tHRD) strategies proposed by Schirris et al.²³. The binary approach involves partitioning the HRD score at the median to distinguish HRP from HRD. That is,

$$mHR{D}^{m}=\left\{\begin{array}{ll}1,&\,\text{if}\,HR{D}^{m}\ge q(0.5)\\ 0,&\,\text{otherwise}\,\end{array}\right.$$

(1)

where HRD^m is the HRD score of the patient m and q(α) indicates the α-quantile of the HRD scores. Another approach utilizes ternary strategy, which divides the HRD scores into three portions. Patients with HRD scores in the top third are assigned HRD, those with HRD scores in the bottom third are assigned HRP, and patients with HRD scores in the middle third are discarded. That is:

$$tHR{D}^{m}=\left\{\begin{array}{ll}1,&\,\text{if}\,HR{D}^{m}\ge q(0.66)\\ 0,&\,\text{if}\,HR{D}^{m}\le q(0.33)\\ \,{\rm{discard}}\, {\rm{patient}}\,,&\,\text{otherwise}\,\end{array}\right.$$

(2)

Compared methods and evaluation metrics

1) Compared methods. For a direct comparison of the performance of SuRe-Transformer with other published methods, we select a weakly supervised method based on Transformer³⁸ and two MIL methods: DMIL³⁹ and DeepSMILE²³. DeepSMILE uses a two-stage MIL, utilizing SSL to pre-train CNN for feature extraction and VarMIL for feature aggregation to achieve HRD status prediction. Only the untrained version of the DeepSMILE code is released to the public. Benchmark testing of these three methods, in comparison with SuRe-Transformer, is conducted on BC pathological images to predict HRD status. All source codes for these methods are downloaded and retrained to ensure a fair and accurate comparison.

2) Evaluation metrics. To evaluate the performance of our framework, we adopt the area under the receiver operating characteristic curve (AUROC, AUC) and F1 score as the primary evaluation metrics, consistent with the previous methods²³. A higher F1 score indicates that the model exhibits high precision and recall in distinguishing between HRD and HRP. This means that the model not only accurately identifies HRD but also minimizes the misclassification of HRP. Additionally, two more curves are adopted for performance analysis, i.e., receiver operating characteristic (ROC) curves with the AUC values and precision recall curves (PRC) with the AUC values.

Implementation detail

We conduct all experiments utilizing a 5-fold cross-validation, where both validation and test sets are split off from the full dataset at the patient level. The dataset is randomly partitioned and is stratified by class, into a training set (80% of cases), and an internal validation (20% of cases). The partitioning was performed at the patient level and therefore all slides from the same patient are always placed into the same set. DINO uses Vision Transformer (ViT) as the backbone and employs the Adam optimizer with a learning rate of 0.0001. For K-means clustering, we take K = 50 following Yan et al.^40,41. Considering computational efficiency, we divide the WSI into four radius ranges (Rad = 0 ~ r, r ~ 2r, 2r ~ 3r, . . . , (m − 1)r ~ mr) in RDSA of SuRe-Transformer, with m = 4. For training SuRe-Transformer, we perform data augmentation. Each WSI is clustered into 50 clusters, and subsequently, patches are then sampled from each cluster based on cluster-size-weighted sampling, resulting in the creation of M patches. Specifically, M takes values of 600, 800, 1000, 1200, 1400, and 1600. Repeating this random selection process S times yields a total of S groups from each WSI, effectively increasing the overall sample size by S times. Here, S = 30. The SuRe-Transformer is trained with AdamW⁴² optimizer with a weight decay of 2 × 10⁻⁵ and a learning rate of 1 × 10⁻⁴. All models are trained for 25 epochs with a batch size of 1. All experiments are performed on a 4-node cluster using the PyTorch framework and scikit-learn. Each node is equipped with a 32 cores CPU, four MI50 GPUs (AMD Vega20, 3840 stream processors, 16GB HBM2 memory), and 128GB of memory.

Overall prediction performance

Given the limited number of studies on HRD status prediction based on WSIs and the variability in labeling methods across different studies, achieving a perfect comparison is challenging. In this paper, we compare the performance of the proposed framework, SuRe-Transformer, with DeepSMILE, DMIL and Transformer-based methods for mHRD and tHRD prediction. To comprehensively demonstrate the predictive performance of our proposed framework, we evaluate the AUROC and F1 scores utilizing a 5-fold cross-validation approach and present the results as mean ± standard deviation. The performance of our proposed framework is shown in Table 1. Our proposed framework with RDSA trained on the patches selected using cluster-size-weighted sampling from WSI achieves the best performance, reaching 0.810 ± 0.047 AUROC and 0.737 ± 0.041 F1 for mHRD, and 0.887 ± 0.034 AUROC and 0.820 ± 0.034 F1 for tHRD. Compared to existing methods, SuRe-Transformer achieves statistically significant improvements across several metrics: it outperforms both DeepSMILE and Wagner et al. in AUROC and F1-score for tHRD (p < 0.05), while it exhibits a significant F1-score improvement over DeepSMILE for mHRD (p < 0.05). The increase in AUROC on mHRD, though numerically higher, is not statistically significant (p > 0.05). Among WSI-based HRD prediction methods, DeepSMILE demonstrates relatively lower performance, primarily due to its reliance on a feature extractor pre-trained on natural images, limiting its ability to effectively capture complex inter-patch relationships characteristic of histopathological images. In contrast, DMIL and Wagner et al. perform better than DeepSMILE as they utilize feature extractors specifically pre-trained on histopathological images. Our proposed method achieves superior performance due to the use of a feature extractor pre-trained on the largest available breast pathology image dataset, effectively capturing relevant histopathological features and inter-patch relationships. Additionally, we visualize the ROC curves and PRC of the above methods in Figs. 3 and 4 for mHRD and tHRD. To evaluate the generalizability of patch embedding based on DINO in our proposed framework, we conduct extensive comparative experiments across all methods. The results indicate that DeepSMILE and DMIL, using DINO for feature extraction, achieve significantly better performance compared to methods that do not incorporate DINO.

Table 1 Comparison of SuRe-Transformer with state-of-the-art methods using the area under the receiver operating characteristic curve (AUROC) and F1 score

Fig. 3: Receiver operating characteristic (ROC) curves and precision recall curves (PRC) analysis of state-of-the-art methods and SuRe-Transformer for HRD status prediction in breast cancer (mHRD).

Top row, from left to right: ROC curves for a DeepSMILE, b Wanger et al., c DMIL and d SuRe-Transformer. Bottom row, from left to right: PRC for e DeepSMILE, f Wanger et al., g DMIL and h SuRe-Transformer. In the ROC curves, the black solid line represents the mean value of 5-fold cross-validation and the shaded area around the curves represents the standard deviation of 5-fold cross-validation. In the PRC, AUC (area under the curve) represents the area under the PRC.

Fig. 4: Receiver operating characteristic (ROC) curves and precision recall curves (PRC) analysis of state-of-the-art methods and SuRe-Transformer for HRD status prediction in breast cancer (tHRD).

Top row, from left to right: ROC curves for a DeepSMILE, b Wanger et al., c DMIL and d SuRe-Transformer. Bottom row, from left to right: PRC for e DeepSMILE, f Wanger et al., g DMIL and h SuRe-Transformer. In the ROC curves, the black solid line represents the mean value of 5-fold cross-validation and the shaded area around the curves represents the standard deviation of 5-fold cross-validation. In the PRC, AUC (area under the curve) represents the area under the PRC.

In addition, we evaluate the effectiveness of RDSA and patch selection based on K-means clustering method with cluster-size-weighted sampling within our proposed framework, using features extracted by DINO_1K. Specifically, when the K-means is not used, patches are randomly selected from each WSI for model training and prediction. As shown in Table 1, the framework utilizing RDSA performs better than the one without it, reaching 0.796 ± 0.030 AUROC and 0.720 ± 0.060 F1 score for mHRD, and 0.869 ± 0.014 AUROC and 0.803 ± 0.035 F1 score for tHRD. We also observe that the framework with cluster-size-weighted sampling increases AUROC by 0.014 up to 0.883 ± 0.028 and increases F1 score by 0.009 up to 0.812 ± 0.032 for tHRD. For mHRD, it demonstrates a comparable performance enhancement, with the F1 score improving by 0.014 to 0.734 ± 0.049. Nevertheless, the improvement in AUROC for mHRD is relatively marginal, exhibiting merely a 0.004 increment to 0.800 ± 0.035. The results illustrate that all components of our proposed framework are essential for the performance of SuRe-Transformer.

Next, we aim to investigate whether morphological features associated with HRD in one cancer type could be used to predict HRD status in other cancer types. We conduct a comparative analysis of the performance of SuRe-Transformer and DeepSMILE by training both models on TCGA-BRCA and subsequently evaluating their effectiveness on TCGA-OV and TCGA-PAAD. The performance of the models is assessed based on AUROC and F1 score for both mHRD and tHRD, as summarized in Table 2. The results show that SuRe-Transformer, trained on TCGA-BRCA, outperforms DeepSMILE in predicting HRD status for TCGA-OV and TCGA-PAAD. Specifically, SuRe-Transformer achieves higher AUROC and F1 scores across all evaluations. For TCGA-OV, SuRe-Transformer achieves an AUROC of 0.713 and an F1 score of 0.664 for mHRD, compared to DeepSMILE’s AUROC of 0.640 and F1 score of 0.471. Similarly, for TCGA-PAAD, SuRe-Transformer achieves an AUROC of 0.599 and an F1 score of 0.577 for mHRD, while DeepSMILE achieves an AUROC of 0.552 and an F1 score of 0.419. Additionally, for tHRD, SuRe-Transformer demonstrates comparable performance improvements over DeepSMILE. These findings indicate that SuRe-Transformer, trained on TCGA-BRCA, can predict HRD status from histological features in other cancer types, suggesting a common “HRD morphology" shared across different cancers. The TCGA-PAAD and TCGA-OV used in this study are available for download at https://portal.gdc.cancer.gov/.

Table 2 Performance comparison of different methods for predicting HRD status on external validation datasets for mHRD and tHRD

Finally, to assess the generalizability and robustness of our proposed framework, we conduct external validation experiments using a total of 192 BC patients (WSI + HRD status) collected from the Beijing ChosenMed Clinical Laboratory Co. Ltd as an independent test cohort. The HRD score is calculated using the scarHRD (https://github.com/sztup/scarHRD), following the methodology outlined in previous studies^43,44. scarHRD uses whole genome sequencing single nucleotide polymorphism (SNP) array data to obtain the patient-wise HRD scores. Based on the 194 WSIs, we conduct an external validation (train on TCGA, test on ChosenMed) and show the results in Table 2. SuRe-Transformer outperforms DeepSMILE in predicting HRD status. Specifically, SuRe-Transformer achieves an AUROC of 0.712 and an F1 score of 0.656 for mHRD, and an AUROC of 0.750 and an F1 score of 0.692 for tHRD. These results highlight the superior performance of SuRe-Transformer in predicting HRD status, underscoring its generalizability and robustness across different sites. However, the predictive performance of SuRe-Transformer on the external validation set drops significantly, however, it is still slightly better than that of DeepSMILE. This is in line with our expectations, because from the naked eye, there is a large difference in color and appearance between the pathological images of ChosenMed and TCGA.

Prediction of breast cancer gene mutation biomarkers

Gene mutation biomarkers are crucial for accurate diagnosis and precise treatment of cancer. In the treatment of BC, TP53 mutations and ERBB2 overexpression^45,46 are associated with poor prognosis and treatment response, while PIK3CA and BRCA1 mutations are linked to the applicability of targeted therapy drugs^47,48. GATA3 mutations⁴⁹ are relevant to the applicability of immunotherapy, and MAP3K1 mutations⁵⁰ are associated with the regulation of cellular signaling pathways. The number of patients with various gene mutations in TCGA BC is shown in Supplementary Table 8. Notably, TP53, PIK3CA, GATA3, and MAP3K1 are the most frequently mutated genes in TCGA-BRCA⁵¹. We evaluate the prediction performance of SuRe-Transformer, which predicts six gene mutation biomarkers on TCGA-BRCA. The results are shown in Table 3. SuRe-Transformer can predict five gene mutation biomarkers (TP53, GATA3, MAP3K1, ERBB2, BRCA1) with an AUROC above 0.844. However, the performance for predicting PIK3CA gene mutation is not as robust. Additionally, we compute the gene variation results on six gene mutations of TCGA-BRCA, as shown in Supplementary Fig. 1. We find that PIK3CA gene mutations, which exhibit poor prediction performance, have higher mutation frequencies, with missense mutations accounting for a significant proportion. The model tends to overfit on missense mutations, leading to poor learning outcomes for rarer mutations, thus affecting the overall predictive performance.

Table 3 Performance evaluation of SuRe-Transformer on 6 different breast cancer gene mutation biomarkers

Additionally, SuRe-Transformer demonstrates superior prediction performance across multiple gene biomarkers compared to the method proposed by Qu et al.⁵², underscoring its robustness and effectiveness. To evaluate the model, we acquire a total of 967 breast cancer patients’ data (WSI + gene mutation status) from TCGA. As shown in Table 4, compared to the method proposed by Qu et al., SuRe-Transformer achieves higher AUROC for all cancer gene biomarkers, including TP53, PIK3CA, MAP3K1, GATA3, ERBB2, and BRCA1. These results indicate that SuRe-Transformer improves prediction performance and provides a more reliable tool for identifying key gene mutations in BC.

Table 4 Performance comparison of SuRe-Transformer and Qu et al.⁵² on different breast cancer gene mutation biomarkers in TCGA-BRCA

The ROC and PRC curves

To comprehensively demonstrate the performance of our proposed framework on TCGA BC dataset, we visualize the ROC curves and PRC of the proposed architecture (Transformer + DINO_1K, Transformer + DINO_1K + Sparse, Transformer + DINO_1K + Sparse + K-means) for tHRD and mHRD classification in Supplementary Fig. 2 and Supplementary Fig. 3, respectively. It can be shown that the framework (Transformer + DINO_1K + Sparse + K-means), which utilizes RDSA and is trained on patches selected using cluster-size-weighted sampling from WSI, achieves the best performance in both ROC curves and PRC. Notably, the PRC for tHRD is relatively flatter compared to that for mHRD, likely because the tHRD discards patients with an HRD score close to the median, leading to a more distinct separation between HRD and HRP.

The visualization of patch clustering result

We follow the experimental settings outlined in previous work from Yan et al.⁴⁰, utilizing K-means clustering with 50 clusters. In this study, we implement the “DINO + K-means” approach to cluster patches within each WSI. The visualization of patch clustering results for BC using K-means (K = 50) is shown in Supplementary Fig. 4. The clustering results indicate that spatially adjacent patches exhibit similar feature representations. Furthermore, discernible spatial heterogeneity of the tumor is observed. Employing clustering-based methodologies facilitates the selection of a diverse array of phenotypic patches, thereby enabling a nuanced characterization of tumor spatial heterogeneity.

Ablation study

We conduct several ablation studies to assess the benefits of different components in our proposed architecture. All experiments are conducted using a 5-fold cross-validation approach on BC dataset from TCGA. First, we evaluate the impact of the number of patches randomly selected from WSI. We fix the rest of the framework (Transformer + DINO_1K), and only change the number of the patches randomly selected from WSI. The experimental results are shown in Supplementary Fig. 5 (random sample). When the input patch sequence length of the Transformer is 1000, it reflects the rich information contained in the WSI. Based on the experimental results, we observe that when the number of randomly selected patches is less than 1000, the performance of the model improves with the increase in the number of selected patches. However, when the number of randomly selected patches exceeds 1000, the performance of the model tends to decrease as the number of patches increases. This phenomenon can be attributed to the longer input sequence, which exacerbates overfitting.

Second, we evaluate the impact of RDSA in ViT. The experimental results are shown in Supplementary Fig. 5 (random sample vs. random sample with RDSA). RDSA effectively increases the input sequence length in ViT. The model achieves optimal predictive performance when the input sequence length is 1200. Based on the experimental results, we observe that when the number of randomly selected patches is less than 1000, the performance of the model with RDSA in ViT is slightly higher than that of the model without RDSA. However, when the number of randomly selected patches exceeds 1000, differences arise in the performance of the model with RDSA compared to that without it. The model with RDSA shows improved performance, whereas the model without RDSA exhibits a decline in performance. This also indicates that RDSA mechanism we introduce effectively alleviates the overfitting issue caused by excessively long input sequences.

Third, we evaluate the impact of the hyper-parameter k in RDSA, with the results reported in Table 5. We fix the rest of the framework (Transformer + DINO_1K + Sparse), and only change the value of the hyper-parameter k within the range (Rad = 0 ~ r, r ~ 2r, 2r ~ 3r, . . . , (m − 1)r ~ mr) in RDSA to observe its effect on model performance, focusing on metrics such as accuracy, F1 score, and AUROC. From the results reported in the Table 5, we observe that the hyper-parameter k with a value of 4, within the range (Rad = 0 ~ r, r ~ 2r, 2r ~ 3r, . . . , (m − 1)r ~ mr) in RDSA, achieves the best performance. Our framework operates within four ranges, resulting in a total of k × 4 = 16. This finding is consistent with the experimental results reported by Yan et al.⁴⁰, which also identifies a value of 16 as optimal. When k is too small, it inadequately represents the input sequence, resulting in suboptimal performance. When k is larger than 4, the performance of the proposed framework does not improve with the increase of k. When k is larger than 8, the performance of the proposed framework is lower than that of not using RDSA. These findings highlight the importance of selecting an optimal k value to balance computational efficiency and model performance in RDSA.

Table 5 Ablation study on the hyper-parameter k in RDSA

Next, we evaluate the impact of the cluster-size-weighted sampling. The experimental results are shown in Supplementary Fig. 5 (random sample with RDSA vs cluster-size-weighted sampling with RDSA). Unlike random patch selection from WSI, cluster-size-weighted sampling ensures that the selected patches are representative. Additionally, we fix the rest of the framework (Transformer + DINO_1K + Sparse + K-means) and only modify the patch selection method: 1). Random patch selection, 2). Patch selection using cluster-size-weighted sampling within each WSI cluster, and 3). Selection of a fixed number of patches from each WSI cluster. The experimental results are shown in Supplementary Fig. 6 (random sample with RDSA vs. kmeans clustering(pro) with RDSA vs. kmeans clustering(equ) with RDSA). The patches selected by cluster-size-weighted sampling exhibit the best performance among the three methods, with a smaller standard deviation. As demonstrated in the systematic studies presented in Supplementary, both random patch selection and cluster-size-weighted sampling from the WSI share the same mathematical expectation, which can approximate the mathematical expectation of all patches in the WSI. However, selecting a fixed number of patches from each WSI cluster, which exhibits the worst performance among the three methods, does not approximate the expectation of all patches in the WSI. The relevant mathematical derivations are provided in Supplementary. Additionally, the patches selected by cluster size-weighted sampling lead to better model performance and a smaller standard deviation compared to randomly selected patches. This is primarily due to the fact that cluster-size-weighted sampling ensures a dispersed selection of patches from the WSI, while random patch selection allows for concentration. The coverage of different patch selection methods on WSIs is presented in Supplementary Fig. 7. As the number of samples increases, cluster size-weighted sampling can cover nearly all patches on the WSI.

Additionally, we evaluate the impact of the clustering embedding module in SuRe-Transformer, with the results reported in Table 6. In our experiments, we fix the rest of the framework (Transformer + DINO_1K + Sparse + K-means), and only change the type of position embedding in SuRe-Transformer. The results show that the clustering embedding module significantly enhances performance compared to other position embedding strategies. This improvement can be attributed to the module’s ability to accurately capture and represent the relative relationships between different patches within the WSIs.

Table 6 Ablation study on the type of position embedding in SuRe-Transformer

Finally, we conduct a comparative analysis of various SSL-based patch embedding methods on pathological images within our proposed framework. Building on recent insights from Wölflein et al.⁵³, we include two SSL methods in the comparison: DINO and CTransPath. Introduced by Wang et al.⁵⁴, CTransPath is a universal feature extractor for pathological images that combines convolutional layers with the Swin Transformer to effectively capture both local and global contextual information. It employs semantically-relevant contrastive learning (SRCL), a novel SSL technique based on MoCo v3⁵⁵, to enhance feature extraction. For a comprehensive comparison, we also include UNI⁵⁶, a universal self-supervised framework that leverages the DINOv2 (ref. ⁵⁷) backbone network and enhances data diversity through pre-training on the “Mass-100K” pathological images dataset, designed for applications such as tissue and disease subtype classification. Additionally, we introduce four variants of DINO: DINO₃₂, DINO₂₀, DINO_1K, and DINO_8K. The DINO₃₂, DINO₂₀, and DINO_1K variants are trained on eight A800 blocks over a span of 36 days, utilizing datasets representing 32 common cancer types, 20 prevalent cancer types, and BC pathological images from TCGA, respectively. Specifically, DINO_8K is trained on 7723 WSIs of BC collected from publicly available datasets, including TCGA-BRCA and CPTAC-BRCA, as well as a private dataset. In our analysis, we maintain the remaining components of the framework (Transformer + Sparse + K-means) and evaluate six patch embedding methods: CTransPath (768D)⁵⁴, DINO₃₂ (384D), DINO₂₀ (384D), DINO_1K (384D), DINO_8K (384D), and UNI (1024D)⁵⁶. The experimental results, presented in Table 7, indicate that DINO_8K achieves the best performance, with an AUROC of 0.887. The comparatively poor performance of UNI can likely be attributed to its exclusion of public pathological image datasets such as TCGA, CPTAC, and PAIP during pre-training. In contrast, DINO and CTransPath utilize datasets like TCGA, which may lead to performance discrepancies, particularly in AUROC results, as shown in Table 7. Additionally, we visualize the attention heatmap of DINO across 6 attention heads in Supplementary Fig. 8. We observe an interesting phenomenon: while foundational models exhibit strong performance across various cancer types, their effectiveness is fundamentally influenced by the volume and specificity of the data used for training. When the focus shifts to a single cancer type, such as BC, it becomes evident that a larger, specialized pre-trained model tailored specifically for BC may yield superior outcomes.

Table 7 The comparison of feature representation learning abilities of different self-supervised learning methods on pathological images

We compare the performance of the proposed framework, SuRe-Transformer, with Transformer-based methods for mHRD and tHRD prediction, including TransMIL, Mamba, and the method proposed by Wagner et al. The performance of our proposed framework is shown in Table 8. SuRe-Transformer achieves the best performance, reaching 0.810 ± 0.047 AUROC and 0.737 ± 0.041 F1 for mHRD, and 0.887 ± 0.034 AUROC and 0.820 ± 0.034 F1 for tHRD.

Table 8 Performance comparison of SuRe-Transformer and Transformer-based method

Discussion

The deployment of HRD biomarker testing in precision oncology for BC presents a promising avenue to enhance treatment efficacy and potentially extend patient survival. However, its adoption has been hindered by significant challenges, primarily due to the complexity and high costs associated with HRD testing, which requires sophisticated diagnostic equipment and specialized expertise. In contrast, tissue sections stained with H&E are ubiquitously available, presenting a viable alternative for assessing HRD status. Therefore, we decide to predict HRD status from WSIs stained with H&E to make HRD testing more accessible and feasible.

In this work, we propose a novel Transformer-based method for the HRD status prediction from WSI, considering the comprehensive correlation among different patches within the WSI. As shown in Fig. 5, there are tens of thousands of patches cropped from a WSI, resulting in extremely long input sequences. The vanilla Transformer exhibits high computational complexity, which restricts the length of the input sequence. Furthermore, the training data for HRD status prediction is limited, making the model susceptible to overfitting. Therefore, using all patches P from each WSI as input is ineffective for predicting HRD status. Instead, we select a small number of informative patches ${P}^{{\prime} }$ and then aggregate them for HRD status prediction. The main contributions are as follows (see Fig. 5):

• To ensure that the selected patches ${P}^{{\prime} }=\{{P}_{k}^{{\prime} }| k=1,\ldots ,S\}$ are informative to preserve the correlation with biomarker labels, we introduce SuRe-Transformer to capture an informative set of patches from WSI. (i) Departing from randomly selecting patches from WSI, we employ cluster-size-weighted sampling to ensure the representativeness of the selected patches. Additionally, we theoretically demonstrate the performance of the framework when using different patch selection methods. (ii) We propose a RDSA mechanism to extend the length of the input sequence in ViT to ensure that sufficient patches are considered. (iii) We pretrain DINO using the largest BC WSI dataset, which consists of 7723 BC WSIs collected from publicly available and private datasets. The resulting model, DINO_8K, serves as an unsupervised feature extractor to obtain better feature representation from image patches within WSIs for clustering and Transformer. (iv) Since the absolute and relative positions between patches after patch sampling are unclear, we propose a clustering embedding module to replace traditional 1D and 2D position embeddings. This module rearranges patches within the same cluster and combines all clusters into SuRe-Transformer’s input sequence, ensuring that spatially or semantically related patches are processed in close proximity.

• Our method achieves better prediction performance. Experiments conducted on BC from TCGA dataset demonstrate the superior performance of SuRe-Transformer in predicting HRD status, compared to the state-of-the-art methods. Additionally, SuRe-Transformer exhibits strong generalizability, effectively predicting HRD status not only in external BC cohorts but also across different cancer types, including ovarian and pancreatic cancers.

• Our framework demonstrates strong generalizability. SuRe-Transformer achieves state-of-the-art performance in predicting key gene mutation biomarkers from BC WSIs. Notably, our proposed framework achieves an AUROC above 0.84 for predicting five key gene mutation biomarkers (TP53, GATA3, MAP3K1, ERBB2, BRCA1).

Fig. 5: Main contributions and key ideas of this article.

Here, ${P}^{{\prime} }=\{{P}_{k}^{{\prime} }| k=1,\ldots ,S\}$ represents the selected patches for analysis in SuRe-Transformer, and S denotes the number of repeated samplings from each WSI.

From a clinical perspective, our method achieves high predictive performance for HRD status based on histopathological slides and demonstrates generalizability across diverse cancer types and patient cohorts, providing a cost-effective and scalable solution. Additionally, it shows promising utility in predicting gene mutations, indicating potential as a versatile biomarker tool, especially in resource-limited settings. Building upon these features, our proposed approach may facilitate the clinical translation of deep learning-based methods.

Although our research has demonstrated promising results, it is essential to acknowledge that our experimental setup has several limitations. First, while SuRe-Transformer demonstrates better generalizability to external patient cohorts compared to DeepSMILE, its performance declines on external test datasets, particularly when compared to internal test datasets. Future work will explore the integration of domain adaptation techniques to address this issue. Second, SuRe-Transformer demonstrates limited effectiveness in handling class imbalance, particularly for underrepresented mutation types. Differences in mutation frequency distributions among various BC gene mutation biomarkers contribute to this imbalance. Although data augmentation techniques were employed to mitigate the issue, the model still struggles with detecting rare gene mutations, especially those with low variant allele frequencies. These limitations may affect its generalizability to larger or more heterogeneous datasets and should be addressed in future work. Third, our current approach has primarily been evaluated using standardized, high-quality datasets such as TCGA, without specific experimentation on lower-quality imaging data. Future studies should investigate the method’s applicability and robustness with lower-quality imaging data to enhance its real-world practicality and effectiveness. Fourth, our proposed method is currently not fully end-to-end, requiring separate preprocessing and feature extraction steps that could introduce variability. Future research will aim to integrate these steps to streamline and enhance the overall workflow. Fifth, the predictive performance of our method exhibits considerable fluctuations. To improve prediction stability and reliability, we intend to introduce a voting mechanism or ensemble approach in future investigations. Finally, the diversity of HRD labeling strategies in existing studies hinders equitable comparisons across different approaches. To provide a reference, we investigate the performance comparisons of existing HRD classification methods based on BC, as presented in Table 9. The HRD status prediction methods utilize various label acquisition methods for HRD status based on the TCGA BC dataset, including HRD score, mutation signatures, and gene signatures. On TCGA BC dataset, SuRe-Transformer achieves similar performance on mHRD classification as Valieris et al.⁵⁸, although a perfect comparison is not available due to the different labeling strategies. Additionally, we observe that the predictive performance of tHRD classification is significantly better than the best reported to date. Specifically, SuRe-Transformer achieves an AUROC of 0.89, while DeepSMILE achieved an AUROC of 0.81.

Table 9 The performance comparison of various methods for predicting HRD status based on TCGA breast cancer

In conclusion, we proposed a novel Transformer-based framework, SuRe-Transformer, for predicting HRD status on WSI from BC with an AUROC of 0.89 ± 0.04 for tHRD and 0.81 ± 0.05 for mHRD. Our model generalizes effectively to an external BC cohort. By publishing the trained models, we hope our proposed framework will provide the researchers and clinicians with new way in HRD prediction using routine histology images.

Methods

Our proposed Transformer-based framework for predicting HRD status comprises four main components: patch embedding based on self-supervised learning, patch selection based on cluster-size-weighted sampling, the architecture of SuRe-Transformer, and radial decay self-attention. Each component will be described in detail in the subsequent subsections.

Patch embedding based on self-supervised learning

In the pre-processing pipeline for WSIs, we use RGB thresholding and Canny edge detection to segment tissue regions, effectively isolating these regions from the background and generating segmentation masks that contain enough tissue cells. Each WSI is tessellated into patches of 256 × 256 pixels (without overlap) at 20 × magnification with a resolution of 0.5 microns per pixel. We view each WSI (known as bag) as a collection of many smaller regions or patches (known as instances). Background patches are removed and the remaining tissue patches are subsequently transformed into feature embeddings. Given the difficulty in obtaining patch-level biomarker annotations, supervised learning methods are not feasible for patch embedding. However, recent advances in SSL offer promising solutions by enabling the training of robust feature extractors without relying on labeled data. Notably, several SSL models trained on large-scale pathology datasets have emerged, providing superior representations for various downstream tasks.

DINO⁵⁹ is a type of SSL method commonly referred to as knowledge distillation, primarily aimed at training smaller networks to replicate the outputs of larger networks for model compression. A comparison of the feature representation learning abilities of different SSL methods is summarized in Table 7. DINO comprises two networks: the student network and the teacher network. The student network is defined by a set of parameters θ_s, as shown in Supplementary Fig. 9. The teacher network has the same network structure as the student network but a different set of parameters θ_t. The teacher network provides the training target for the student network, whose parameters are updated by using an exponential moving average (EMA) on the student weights, following the rule:

$${\theta }_{t}\leftarrow \lambda {\theta }_{t}+(1-\lambda ){\theta }_{s}$$

(3)

where λ follows a cosine scheduled from 0.996 to 1 during training. Given a fixed teacher network g_θt, we learn to match the distributions by minimizing the cross-entropy loss with respect to the parameters of the student network g_θs,

$$\mathop{\rm{min}}\limits_ {{\theta }_{s}}\,H({P}_{1},{P}_{2})$$

(4)

where $H(a,b)=-a\log b$. After training converges, the patch embedding output g_θt is used as the pathological image representation. Therefore, we utilize DINO, which has been pre-trained on 6 cohort of BC pathological images, as shown in Fig. 1a, to convert each image patch of size 256 × 256 × 3 into a discriminative feature vector of size 384. We save each patch in a hierarchical data format hdf5, which contains an array of extracted features along with their corresponding patch coordinates.

Patch selection based on cluster-size-weighted sampling

Transformer-based architectures are computationally expensive due to the self-attention mechanism, the computation of which is typically quadratic to the length of the input sequence. This poses a considerable challenge for analyzing histopathological images, such as WSIs, which are of extremely high resolution, thereby limiting the number of input patches in existing Transformer-based WSI analysis approaches. Traditional Transformers are only capable of handling short sequences of less than 1000²⁸. To address this limitation, we select several patches from WSI and aggregate these patches into WSI-based feature representations. However, random patch selection is insufficient for capturing the representative information from the WSI.

Ensuring the representativeness of the selected patches is essential for accurate HRD status prediction. Therefore, we propose a patch selection method based on cluster-size-weighted sampling, as shown in Fig. 1b, where all patches cropped from the WSI are clustered into K clusters using the K-Means algorithm, a standard unsupervised clustering technique. Subsequently, patches are sampled from each cluster based on the proportion of patches in each cluster relative to the total number of patches in the WSI, ensuring that the required number of patches, denoted as M, is obtained from each WSI. Detailed information can be found in Section “SuRe-Transformer architecture”. Specifically, each patch is transformed into a descriptive, one-dimensional feature representation with a length of 384 based on the DINO_8K. The K-means algorithm is then applied to cluster these feature vectors into K clusters. Patches are sampled from each cluster in accordance with the aforementioned selection proportion, ensuring that the selected patches effectively represent the diversity within each cluster for subsequent analysis.

SuRe-Transformer architecture

There are two main obstacles that hinder the application of Transformers for HRD status prediction. First, the inherent computational complexity of the vanilla Transformer imposes limitations on the length of the input sequence. Second, the limited availability of datasets for HRD status prediction increases the risk of model overfitting. Therefore, it is crucial to explore alternative Transformer architectures tailored to address these challenges.

The SuRe-Transformer builds upon the ViT. We define the i^th patch obtained from gigapixel WSIs as an instance x_i, and consider the set of all patches in the WSI as a bag (B = {x_i∣i = 1, . . . , N}), where N represents the number of patches in the WSI. The input of SuRe-Transformer architecture is a dataset P consisting of a sequence of M patches, with the dataset label derived from the corresponding WSI label:

$$P=\{{P}_{k}^{{c}_{k}}| k=1...K,{c}_{k}=1...{C}_{K}\}$$

(5)

where K is the number of clusters per WSI, ${C}_{k}=\frac{{N}_{k}}{N}\times M$ denotes the number of selected patches from the k^th cluster, M satisfies the following relationship $M=\mathop{\sum }\nolimits_{k=1}^{K}{C}_{k}$, and N_k is the number of patches in the k^th cluster. The patch embedding E(⋅) is extracted based on the DINO_8K.

Since the Transformer does not account for the position embedding³⁸, position embeddings are added to the patch embeddings to capture the order of the input sequence. However, given that the input of the SuRe-Transformer is sparsely distributed, traditional 1D and 2D position embeddings fail to effectively capture the relative relationship between the patches within the WSI. To address this issue, we propose a clustering embedding module in SuRe-Transformer, making it more suitable for the HRD status prediction. As illustrated in Supplementary Fig. 10, the vanilla Transformer encoder utilizes [*, Patch embedding] plus position embedding. In contrast, the SuRe-Transformer encoder introduces a clustering embedding module, which arranges the patches of the same cluster together and combines all clusters to generate the input sequence of SuRe-Transformer, ensuring that spatially or semantically related patches are processed in closer proximity. This allows the SuRe-Transformer to better understand the contextual information, as each patch is assigned its embedding information based on its position in the sequence.

As illustrated in Fig. 1c, each layer of the SuRe-Transformer encoder consists of a multi-head radial decay self-attention (MHRDSA) module, which is detailed in the subsequent subsection, and a multi-layer perceptron (MLP) block. We follow the ViT, adding the CLS token at the beginning of the input sequence P. The input token embeddings can be represented as:

$${z}_{0}=[E(CLS);E(P)]$$

(6)

The Transformer encoder module consists of L-stacked layers, where LN denotes the layer normalization, and MLP comprises two fully connected layers with a non-linear GELU activation function²⁵. The SuRe-Transformer can be expressed as follows:

$${z}_{l}^{{\prime} }={MHRDSA}(LN({z}_{l-1}))+{z}_{l-1},l=1,...,L$$

(7)

$${z}_{l}=MLP(LN({z}_{l}^{{\prime} }))+{z}_{l}^{{\prime} },l=1,...,L$$

(8)

The CLS token from the output of the L Transformer layers is used as the WSI representation and is passed into the MLP head for predicting the HRD status for BC. The prediction of the HRD status can be formulated as follows:

$${HRDstatus}=MLP({z}_{L}[0])$$

(9)

Radial decay self-attention

The self-attention (SA) in the Transformer computes the query-key product as follows:

$$Att(Q,K,V)=softmax\left(\frac{Q{K}^{T}}{\sqrt{{d}_{k}}}\right)V$$

(10)

where Q, K, and V represent the queries, keys, and values, respectively, and d_k denotes the dimensions of Q and K. The multi-head self-attention (MHSA) applies SA in every head and concatenates the heads in a weighted manner:

$${MHSA}=concat(hea{d}_{1},...,hea{d}_{h}){W}_{o}$$

(11)

where head_i = Att(Qⁱ, Kⁱ, Vⁱ), i = 1. . . h, and W_o are learnable parameters. Specifically, we choose a Transformer network architecture consisting of six layers, each with six heads. The MHSA is an extension of the SA mechanism in Transformer that enables the model to learn in different representation spaces by parallelizing multiple independent SA operations, allowing it to more effectively capture complex relationships within the input sequence. Most ViT-based architectures, such as the original ViT, perform poorly on small-scale biomarker datasets due to a lack of inductive bias. In the task of HRD status prediction, it is unnecessary for each patch in the input sequence to establish relationships with all other patches. Instead, it suffices for each patch to interact with only a subset of key patches.

We introduce a RDSA mechanism, inspired by Generative Pre-trained Transformer 3 (GPT-3)⁶⁰. One fundamental idea behind GPT-3 is sparse self-attention, where attention is restricted to tokens within a relative distance of no more than k. That is, attention is selectively applied to tokens at relative distances of k, 2k, 3k, and beyond, while remaining zero for all other tokens. While the concept based on GPT-3 is one-dimensional, pathological images are two-dimensional; hence we analogize the distance (K = k, 2k, 3k, . . .) in 1D context to the radius (Rad = 0 ~ r, r ~ 2r, 2r ~ 3r, . . . , (m − 1)r ~ mr) in 2D context. The illustration of the RDSA mechanism is depicted in Supplementary Fig. 11. The key idea is to reduce the dense attention matrix to a sparse representation by restricting attention computation to a limited subset of query-key (QK) pairs.

In this paper, we first obtain the attention matrix based on patches from the model, as the [CLS] token does not contain WSI coordinate-related information. Subsequently, we compute the relative distances D between patches within WSIs, as illustrated below:

$${D}_{ij}=\sqrt{{({x}_{i}^{* }-{x}_{j}^{* })}^{2}+{({y}_{i}^{* }-{y}_{j}^{* })}^{2}}$$

(12)

$${x}_{i}^{* }={x}_{i}/\max (x)$$

(13)

$${y}_{i}^{* }={y}_{i}/\max (y)$$

(14)

where (x_i, y_i) represents the coordinate of the i^th patch in the WSI, and max(x) denotes the maximum value of the coordinates x among all patches. We first mark the patches within each radius range by:

$${D}_{ij}^{\left(n\right)}=\left\{\begin{array}{ll}1,&\,{\text{if}\, (n-1)}\,r < {D}_{ij} < nr\\ 0,&\,\text{otherwise}\,\end{array}\right.,\quad n=1,2,3,\ldots ,m$$

(15)

To implement the patch selection within the radius range, we utilize a randomly initialized unit matrix R within the (0, 1) range to assign a selection probability to the patches within the different radius ranges, as illustrated by S.

$${S}^{(n)}=R\otimes {D}^{(n)},n=1,2,3,...,m$$

(16)

where ⊗ denotes element-wise operation. Let ${t}_{i}^{(n)}$ denote the k^th largest values of range n in the i^th row, and k is the hyper-parameter. To select the top-k values of S in the different radius ranges Rad, we use the mask operation M(.), which is illustrated as:

$${M}_{ij}^{(n)}=\left\{\begin{array}{ll}{P}_{ij},&\,\text{if}\,{S}_{ij}^{(n)}\ge {t}_{i}^{(n)}\\ -\infty ,&\,\text{otherwise}\,\end{array}\right.$$

(17)

We concatenate the {M⁽ⁿ⁾; n = 1, 2, …, m} to form M = [M_ij]:

$${M}_{ij}=\left\{\begin{array}{ll}{M}_{ij}^{(n)},&\,\text{if}\,\exists {M}_{ij}^{(n)}\,\ne -\infty \\ -\infty ,&\,\text{if}\,\forall {M}_{ij}^{(n)}=-\infty \end{array}\right.$$

(18)

After concatenation, the RDSA is defined as follows:

$${RDSA}=softmax(M)V$$

(19)