An interpretable deep learning model for detecting BRCA pathogenic variants of breast cancer from hematoxylin and eosin-stained pathological images

Patient selection

Medical data from breast cancer patients who underwent BRCA gene testing at Chongqing University Cancer Hospital (CUCH) between May 2017 and December 2022 were retrospectively reviewed. The inclusion criteria were: (1) Diagnosis with primary breast invasive ductal cancer (IDC) or invasive Lobular carcinoma (ILC); (2) conduct of BRCA gene or genomic testing on blood and/or tumor samples, which included the examination of germline and/or somatic variations; (3) availability of operative pathological H&E images before antineoplastic treatments were available; The exclusion criteria were: (1) Lack of operative pathological H&E images; (2) poor quality of H&E-stained images; (3) presence of bilateral, multifocal, or special invasive breast cancer. Adhere to the WHO guidance on Ethics & Governance of Artificial Intelligence for Health was maintained. The study was approved by the Ethics Committee of Chongqing University Cancer Hospital (Ethics number: CZLS2023213-A) and patient consent was waived for this retrospective analysis. In reporting the results, the Standards for Reporting Diagnostic Accuracy Studies (STARD) 2015 guidelines were strictly followed.

This study included two cohorts. The first cohort consisted of 152 patients (217 WSIs) collected from Chongqing University Cancer Hospital between May 2017 and December 2022. The second cohort consisted of 102 patients (102 WSIs) collected from The Cancer Genome Atlas (https://gdc.cancer.gov/). The CUCH cohort was randomly divided into a training set and an internal test set, while the TCGA cohort served as an independent external test set. The training set was utilized for hyperparameter tuning through cross-validation, and the test sets were used to assess generalization performance.

Genomic DNA samples were extracted from peripheral blood and/or surgical tissue samples. The BRCA genetic testing was conducted using next-generation sequencing (NGS) technology. Variants were annotated following the Human Genome Variation Society (HGVS) nomenclature guidelines (http://varnomen.hgvs.org/). The biological significance of all reported variants was assessed using the ClinVar database (www.ncbi.nlm.nih.gov/clinvar/). According to the guidelines of the American College of Medical Genetics and Genomics (ACMG) (Richards et al., 2015), the detected variants were classified as pathogenic variants, likely pathogenic variants, variants of uncertain significance, likely benign, or benign. Patients identified with pathogenic and likely pathogenic variants were categorized into the BRCA PV group, which includes BRCA1 PV and BRCA2 PV. Patients identified with variants of uncertain significance, likely benign, or benign variants were categorized into the BRCA WT group. Pathogenicity was defined as any copy loss event (heterozygous or homozygous) in the BRCA1 and BRCA2 genes, as well as truncating or exonic single-nucleotide alterations/indels less than 10 bp in length that are labeled as “pathogenic/likely pathogenic” in ClinVar. The study focused on predicting these two binary outcomes, BRCA PV and BRCA WT for breast cancer patients.

Image preprocessing and sample preparation

All samples were fixed in 4% neutral formalin, embedded in paraffin, cut into 4 μm thick sections, and stained with H&E. The H&E-stained histopathological slides were scanned at 40× magnification using a KFBIO KF-PRO-005 digital scanner and saved in SVS format. Experienced pathologists with over 1 year of specialization in breast cancer pathology, performed quality control on the images. Slides of poor quality, particularly those with tissue folds and blurriness, were excluded. The TCGA cohort exclusively utilized 40× FFPE WSIs labeled as diagnostic, while excluding slides marred by pen marks or inadequate staining. The WSIs were segmented into tiles of 512 × 512 pixels, with 50% overlap. Tiles displaying more than 75% of the background were discarded. Observed color variations, attributable to differences in raw materials, manufacturing processes, staining methods, and different digital scanners, can cause models to focus on color differences rather than essential tissue morphology for analysis. To address this issue, color normalization was applied to all tiles using the structure-preserving color normalization (SPCN) technique (Vahadane et al., 2016). The model standardized the color pattern of all tiles to resemble that of the target (Fig. S1). Various data augmentation techniques were employed, including random flipping, random rotation, cropping, and adjustments to brightness, contrast, saturation, and hue.

A stratified sampling approach was used to randomly assign images to the training and internal test set, maintaining an 8:2 ratio at the patient level. This method ensures that images of the same patient are assigned exclusively to either the training set or the internal test set, preventing images from the same patient from appearing in both sets simultaneously.

Segmentation network

A segmentation network was constructed to automatically identify tumor regions in WSIs. Two experienced pathologists, each with over 3 years of experience in breast cancer pathology, annotated tumor areas on 145 WSIs using the open-source software Qupath Image Scope. These annotated areas were subsequently reviewed independently by a pathologist with over 10 years of experience. The labeled slides were randomly divided into a training set (80%) and a test set (20%), ensuring that WSIs from the same patient did not appear in both sets. Utilizing the ResNet 34 architecture, the network processed tiles from both tumor and non-tumor regions. Deep convolutional and pooling layers extracted features from each tile, which were then passed to the fully connected layer. Following the fully connected layer, a softmax function generated a probability (tumor or non-tumor) for each tile. The network’s performance was evaluated by comparing its results with the manual annotations made by three pathologists of varying experience levels-junior, medium, and senior, with 2, 5 to 8, and 10 years of experience, respectively. Ultimately, this segmentation network was utilized to automatically segment tumor regions in the TCGA cohort.

Prediction network

The CUCH cohort comprised 217 WSIs from 152 patients, each contributing one or two slides. The dataset was stratified and randomly divided into a training set (633,484 tiles from 174 WSIs) and an internal test set (112,096 tiles from 43 WSIs). The BiAMIL network was developed using a multi-instance learning (MIL) framework and attention mechanisms. This architecture included three main components: the feature extraction module based on the ResNet 34 model (He et al., 2016), the bi-directional self-attention module, and the classification module. A bag containing N image tiles was inputted into the feature extractor, generating a feature matrix with dimensions of N × 1000. This feature matrix was processed through the attention module for aggregation. In the bi-directional self-attention module, data features were transformed into embedding vectors via a feature embedding layer. Two attention heads were then designed: a high-risk PV head and a low-risk PV head. Each head utilizes independent core weight matrices $W^k,W^q,W^v$. For each head, the 1,000-dimensional embedding features in each instance are transformed into matrices $K=\left\{W^k{i}_1,W^k{i}_2,\cdots ,W^k{i}_N\right\}$, $Q=\left\{W^q{i}_1,W^q{i}_2,\cdots ,W^q{i}_N\right\}$ and $V=\left\{W^v{i}_1,W^v{i}_2,\cdots ,W^v{i}_N\right\}$ through linear transformations. The attention score matrix $\alpha$ is defined as $\alpha =K^{T}Q$. Each element ${\alpha }_{i,j}$ in $\alpha$ is computed as the inner product of the i-th row of $Q$ and the j-th column of $K^T$, for example ${\alpha }_{1,2}=k_1\cdot q_2={\sum }_{i=1}^D{k}_{1,i}{q}_{2,i}$. Subsequently, the attention matrix A is calculated as $A=V\times Softmax\left(\alpha \right)$. For each embedding feature vector, the corresponding attention vector is denoted by $a_i=A_{v_i}$, where $v_i$ represents the i-th column of matrix $V$, such that $a_i={\sum }_{n=1}^N{v}_i{\alpha }_{i,n}$. Let $A_h$ and $A_l\in R^{N\times 1000}$ denote the attention matrices generated separately by the two-head self-attention mechanism. The overall attention matrix $A_c={\left\{A_h:A_l\right\}}_h$ is formed by concatenating $A_h$ and $A_l$. A weight vector w is used to transform $A_c$ into the final attention weights $F_{AW}^T\in R^{1000\times 1}$ via a linear transformation $F_{AW}=\sigma \left(W_A^T{A}_c\right)$. The resulting $F_{AW}^T$ represents the attentional representation of the instance bag. The weight vector $W_A$ plays a critical role in governing the learning process of the high-risk and low-risk PV heads, optimized using the following loss function:

$$loss=-\sum _{i=1}^C{y}_{i}log\left(Softma{x}_{MLP}\left(\gamma mask{\left(W_A^T{A}_c\right)}_h+\mu mask{\left(W_A^T{A}_c\right)}_l\right)\right).$$

The hyperparameters $\gamma$ and $\mu$ are initially initialized to 0.5. As the learning process progresses, these parameters gradually converge towards their optimal values. Two sets of attention-weighted discriminative feature vectors were generated, each set having a dimensionality of 1,000. These two 1,000-dimensional integrated feature vectors were merged and inputted into the classification module. The classification module consisted of a three-layer multilayer perceptron (MLP) with 1,000, 64, and two neurons in each layer, respectively, utilizing the Tanh function as the activation function. The prediction probability for each instance bag was ultimately output through the softmax function. To optimize the model, we employed a five-fold cross-validation by randomly dividing the training set into five balanced subsets (Fig. S2). Using the Adam optimizer, the optimal hyperparameters were identified as a bag of N = 35, a learning rate of 0.001, a batch size of 18, and total epochs of 20 (with the learning rate decreasing by 0.1 every ten epochs). At last, the network was evaluated on the internal and external test sets.

Tissue features visualization and cell features quantitative analysis

To elucidate the critical histological features most influential in the BiAMIL model’s predictions, the top 20% of tiles, based on their attention scores for predicting BRCA status, were selected. Attention-based visualization techniques were employed to visually represent these tiles. First, the weight values for each tile were calculated. Subsequently, the top 20% of tiles in WSIs with high attention weights were identified, and Smooth Grad-CAM was used to pinpoint the regions within each tile that the neural network relied on for generating predictions. We used cell detection and classification functions using QuPath to identify three types of cells, tumor cells, tumor-infiltrating lymphocytes (TILs), and stroma. TILs% calculated as follows:

$$TILs\mathrm{\%}={\displaystyle \frac{TILs}{TILs+Tumorcells}\times 100\mathrm{\%}.}$$

Cell features, including cell morphology, nuclear morphology, and staining features, were extracted from the tumor regions (Bankhead et al., 2017). Cells were then classified into different clusters using the k-means clustering method based on their cell features. These clusters were projected onto a 2-D UMAP projection (Dorrity et al., 2020). The differences in cell features between BRCA PV and WT were assessed using the Wilcoxon rank-sum test.

Statistical analysis

The performance of the model was evaluated at the slide level. A probability value was generated for each bag, and a final probability for each slide was determined by calculating the average probability across all bags within the slide. The model’s performance was assessed using several measures, including area under the receiver operating characteristic curve (AUROC), accuracy, sensitivity, specificity, positive predictive value (PPV), negative predictive value (NPV), and F1 score. The most relevant cell features for discriminating between BRCA PV and BRCA WT were identified using the Wilcoxon rank-sum test. All statistical analyses were performed using scikit-learn 0.24.2 in Python and R (version 4.1.1; R Core Team, 2021), with a p-value of less than 0.05 was considered statistically significant.

Clinical characteristics of patients

A total of 254 patients from two independent cohorts were included in our study based on the inclusion and exclusion criteria. Specifically, 152 patients from the CUCH cohort were randomly assigned to a training set and an internal test set, while 102 patients from the TCGA cohort served as an external test set. Detailed clinicopathological characteristics of the patients across the training, internal test, and external test sets are presented in Table 1. In the training set, 62 (50.8%) of patients were ER and PR positive, 12 (9.8%) were HER2 positive, and 48 (39.4%) were triple-negative. In the internal test set, 16 (53.4%) of patients were ER and PR positive, five (16.6%) were HER2 positive, and nine (30.0%) were triple-negative. In the external testing set, 70 (68.6%) of patients were ER and PR positive, 11 (10.8%) were HER2 positive, and 18 (17.6%) were triple-negative. In the training set, 29 patients had BRCA PV, with 17 (58.6%) for BRCA1 and 12 (41.3%) for BRCA2. The internal test set included seven patients with BRCA PV, comprising 4 (57.1%) with BRCA1 and three (42.9%) with BRCA2. In the external test set, among 26 patients with BRCA PV, 13 (50.0%) had BRCA1, 11 (42.3%) had BRCA2, and two (7.7%) had both.

Performance of the segmentation network

A segmentation network was developed to accurately distinguish tumor regions from WSIs. This model demonstrated high efficacy, achieving an AUC of 0.960 (95% CI [0.959–0.961]). The accuracy, sensitivity, and specificity were 0.888 (95% CI [0.887–0.890]), 0.859 (95% CI [0.856–0.861]), and 0.908 (95% CI [0.906–0.909]), respectively. Compared with three pathologists of varying experience (junior, medium, and senior), the model’s performance was almost equivalent to that of a medium-level pathologist (Fig. 2A). The confusion matrix shows specific performance metrics of the segmentation model. Among tumor region samples, the model correctly identified 86% as tumor, while misclassifying 14% as non-tumor. Among non-tumor region samples, 91% were correctly identified as non-tumor, with 9% misclassified as tumor (Fig. 2B). Three original WSIs and their corresponding annotation maps are illustrated in Fig. 2C, where the red areas in the maps denote the tumor regions identified by the model. The alignment of the model-generated annotation maps with the tumor regions in the original images validates the model’s accuracy.

Performance of BiAMIL network

The BiAMIL model was developed by first inputting normalized tiles into the ResNet 34 model to extract 1,000-dimensional features. These features were then aggregated using a Bi-directional self-attention mechanism approach. Subsequently, the aggregated feature matrix was fed into a multilayer perceptron (MLP) classifier to predict the probability of BRCA gene mutation status for each WSI (Fig. 3). BiAMIL achieved AUC values of 0.819 (95% CI [0.673–0.965]) in the internal test set (Fig. 4A), and 0.817 (95% CI [0.712–0.923]) in the external test set (Fig. 4B). The five-fold cross-validation results for the BiAMIL model on the training set are displayed in Table S1. The model was compared with the ResNet 34 by Wang et al. (2021). ResNet 34 achieved AUC values of 0.783 (95% CI [0.624–0.941]) the internal test set (Fig. 4C), and 0.684 (95% CI [0.559–0.810]) in the external test sets (Fig. 4D). The five-fold cross-validation results for the ResNet 34 model on the training set are shown in Table S2. In the internal test and external test sets, BiAMIL outperformed the ResNet 34 model in accuracy, sensitivity, specificity, positive predictive value (PPV), negative predictive value (NPV), and F1 score (Figs. 4E and 4F). Due to hormone receptor-positive and triple-negative breast cancer accounting for 80% of breast cancer cases and the high occurrence of BRCA1 in triple-negative breast cancer, a subgroup analysis of BiAMIL’s performance in these subtypes was conducted, with results shown in Table 2. These findings demonstrate that BiAMIL is a more accurate and reliable method for estimating the BRCA status compared to the traditional ResNet 34 model.

Interpretability analysis of tissue features

To explore the relationship between pathological tissue features and BRCA gene status, a post-hoc explanation of the BiAMIL model was conducted using the Class Activation Mapping (CAM) algorithm. The model generated heatmaps for tiles with high attention scores, indicating that tumor areas, including tumor cells and tumor-infiltrating lymphocytes (TILs), are significant contributors to the model’s predictions, while the stroma contributes less (Figs. 5A and 5B). Tiles predicted as BRCA PV exhibited a denser number of TILs compared to those predicted as BRCA WT, as evidenced by the median percentage of TILs being higher in the BRCA PV group than in the BRCA WT group (25.8% vs. 14.1%) (Fig. 5C). The effectiveness of the BiAMIL model in distinguishing BRCA gene mutation status was further assessed using the t-SNE algorithm, revealing a clear separation between the BRCA PV and BRCA WT groups (Fig. 5D).

Interpretability analysis of cell features

To further explore the relationship between tumor cell features and BRCA gene status, the Qupath software was used to automatically identify tumor cells in regions highly influenced by the BiAML model. Pathologists reviewed and refined the delineation of tumor cells in each WSI, extracting features of cell morphology. The cells were then clustered into six groups using the k-means, with each cluster represented by a different color (Fig. 6A). A correlation analysis was constructed between the clusters and cell morphological features, illustrating the relationship between each cluster (Cluster 0 to Cluster 5) and the cell morphological features. This analysis indicated a strong correlation between Cluster 2 and the morphological features of tumor cells and nuclei (Fig. 6B). Subsequent analysis of the differences in tumor cell and nuclear morphological features between the BRCA PV and WT groups within Cluster 2 revealed that, compare to the BRCA WT group, the BRCA PV group exhibited larger tumor cell areas (501.75 vs. 461.50, P = 1.41e−13), larger nuclear areas (180.00 vs. 168.75, P = 3.13e−09), and higher nuclear cell area ratio (0.36 vs. 0.34, P = 5.73e−11) (Figs. 6C–6E).