Direct prediction of Homologous Recombination Deficiency from routine histology in ten different tumor types with attention-based Multiple Instance Learning: a development and validation study

Data Acquisition

In total data from 5,155 patients of 10 tumor types from The Cancer Genome Atlas (TCGA) and 573 patients from five tumor types from the Clinical Proteomic Tumor Analysis Consortium (CPTAC, Figure 1C) were obtained from https://www.cbioportal.org/. Accordingly, the cancer types included in the present study were breast invasive carcinoma (TCGA-BRCA n=1,058), colorectal cancer (TCGA-CRC n=580), glioblastoma (TCGA-GBM n=420, CPTAC-GBM n=99), liver hepatocellular carcinoma (TCGA-LIHC n=364), lung adenocarcinoma (TCGA-LUAD n=536, CPTAC-LUAD n=111), lung squamous cell carcinoma (TCGA-LUSC n=497; CPTAC-LSCC n=109), ovarian cancer (TCGA-OV n=520), pancreatic adenocarcinoma (TCGA-PAAD n=177; CPTAC-PDA n=153), prostate adenocarcinoma (TCGA-PRAD n=488) and endometrial carcinoma (TCGA-UCEC n=515, CPTAC-UCEC n=101, Supplementary Figure 1A,B). Image data and clinical data were available in TCGA-BRCA for n=1005, TCGA-CRC for n=496, TCGA-GBM for n=232, CPTAC-GBM for n=99, TCGA-LIHC for n=348, TCGA-LUAD for n=460, CPTAC-LUAD for n=106, TCGA-LUSC for n=451, CPTAC-LSCC for n=108, TCGA-OV for n=90, TCGA-PAAD for n=173, CPTAC-PDA for n=139, TCGA-PRAD for n=391, TCGA-UCEC for n=467 and CPTAC-UCEC for n=99, therefore leaving us in total with n=4,565 patients for the analysis (Figure 1C, Supplementary Figure 1A,B). Moreover, some figures were created using https://www.cbioportal.org/ (37,38). For additional experiments on BRCA1/2 mutational status we retrieved data from Riaz et al previously published paper (39). Estrogen receptor data for the subgroup analysis was only available for n=661 patients in the TCGA-BRCA cohort.

Image Preprocessing

WSIs were downloaded for the above mentioned cohorts from the GDC Portal (https://portal.gdc.cancer.gov/) and The Cancer Imaging Archive (TCIA, https://www.cancerimagingarchive.net/). Initially, the images were tessellated into patches with an edge length of 256 μm and a resolution of 224×224 pixels. Secondly, the patches for each cohort were color normalized using the Macenko spectral matching technique(40) to enforce a standardized color distribution across cohorts. To train the prediction models, we used our in-house open-source DL pipeline “marugoto” (https://github.com/KatherLab/marugoto) consisting of a self-supervised learning (SSL) model using a pre-trained ResNet50 architecture with ImageNet weights, fine-tuned pan-cancer on approximately 32.000 WSI to extract a 2048-dimensional feature vector for each patch per patient (41). To obtain patient-level predictions, 512×2048 feature matrices (MIL bags) are constructed by concatenating 512 feature vectors selected at random per patient and fed into an attMIL framework with the following architecture: (512×256), (256×2) with a subsequent attention mechanism (Figure 1B). (42,43)

Calculation of HRD Scores

For the patient-wise calculation of HRD, single nucleotide polymorphism (SNP) data, generated by the Allele-Specific Copy number Analysis of Tumors (ASCAT) algorithm, was downloaded from the Genomic Data Commons (GDC) Portal: https://portal.gdc.cancer.gov/ (accessed 06/15/2022). In CPTAC, the respective data was only available for the CPTAC-3 cohort. The HRD score was calculated using the scarHRD (https://github.com/sztup/scarHRD), as described in previous studies (20,44). ScarHRD uses whole genome sequencing data in the form of SNP arrays to calculate the three subscores LOH, LST and TAI. The sum of these subscores makes up the patient-wise HRD score (Figure 1A). The cut-offs of the different subscores have been previously defined by Abkevich et al. for LOH, Popova et al. for LST and Birkbak et al. for TAI (16–18). Adding up the LOH, LST and TAI scores, patients can be divided into HRD high (HRD-H) and HRD low (HRD-L) at a cut-off of 42 (7). All patients in the CPTAC-GBM cohort were HRD-L. Hence, we excluded them from further analysis (Supplementary Figure 1A,B).

Experimental Design

In our study, we performed three main experiments (Figure 1B). To assess the baseline predictability of HRD from routine histology, we first trained an HRD classifier in a within-cohort approach using five-fold-cross-validation within each of 10 tumor entities mentioned above in the TCGA cohorts (internal validation). This was achieved by randomly splitting each cohort on the level of patients, creating non-overlapping training and test sets for model training. The ratio for splitting the training and testing set was 80:20 of the entire dataset, and the training and validation set was split 75:25 of the training set. Thus, the absolute split for training, internal validation and internal testing was 60, 20 and 20, respectively. Five different models were trained until each part was used as a test set once. Thus, no data leakage from the test set to the training set occurred. This process was repeated individually for each cancer type in the TCGA cohorts. A weighted cross-entropy loss function was used to assist the model with the imbalanced dataset. Secondly, we deployed the five models trained in the first experiments on the same tumor type from the CPTAC cohorts as an external validation. By utilizing this approach, we circumvent any potential claims of selecting the model with the highest AUROC in the external validation. Lastly, we trained an HRD classifier on the TCGA-BRCA cohort, which had the highest number of patients, and deployed it on all other TCGA cohorts (CRC, GBM, LIHC, LUAD, LUSC, PRAD, PAAD, OV, UCEC) as well as on all CPTAC cohorts (LUAD, LSCC, PDA, UCEC). In our study, we aimed to evaluate the performance of the models using the AUROC, which is commonly used for assessing the accuracy of binary classification tasks. Our primary statistical endpoint was the AUROC +/-95%-confidence interval (CI) and Area under the precision recall curve (Supplementary Table 1). To further assess the performance of each model, we used a two-sided t-test to compare the patient-level prediction scores between the HRD-H and HRD-L patient groups as defined by the ground truth and report the p-values, assuming a significance level of 0.05 as statistically significant, without correction for multiple testing (Supplementary Table 1). As a final step to obtain a more profound understanding of the TCGA-BRCA cohort, we uploaded our custom HRD-H and HRD-L ground truth and predicted subgroups in cbioportal to examine the characteristics of these cases in the TCGA-BRCA PanCancer Atlas cohorts.

Explainability

To visualize the output of our model, we created high resolution heat maps that show the spatial distribution of our model’s attention and prediction scores on the WSI. Therefore, using RetCCL convolutional neural network image feature vectors for 32×32 pixel fields were extracted from the WSI. We then calculated attention and classification scores for each image region and normalized them across the distribution of scores within each patient cohort. Based on these scores, color heatmaps for each patient, with red indicating high attention or a positive classification and blue indicating low attention or a negative classification were generated. To ensure interpretability of the underlying morphology together with the attention and classification scores, we separately reconstructed the final attention and classification heatmaps by blending the raw color heatmaps with the image features. This approach allows us to interpret the output of our model in a way that is easy to understand and provides insight into the underlying morphology of the tumor.

HRD is predictable from histology with attMIL

First, we investigated whether DL could predict HRD status from H&E types within 10 different types of cancer from the TCGA cohort. We used cross-validation on the level of patients to train and test an attMIL-based DL model within each cohort. In our dataset, the prevalence of HRD ranged from 3% in glioblastoma (GBM) up to 63% in OV (Supplements Figure 1C). We found that in five out of 10 cancer types, the mean prediction AUROC was above 0.6, and the 95% CI of the fold-wise HRD prediction AUROCs remained above the null hypothesis of 0.5. Among these, HRD prediction reached statistical significance with a p-value below 0.05 in three cancer types: endometrial cancer (UCEC, AUROC 0.79+/-0.04, p=0.0008), breast cancer (BRCA, AUROC 0.78+/-0.02, p<0.0001) and lung adenocarcinoma(LUAD, AUROC 0.66+/-0.05, p=0.02; Figure 2A). AUPRC values are reported in the Supplementary Table 1. Prediction of HRD was not possible in LUSC, LIHC, GBM, as their prediction AUROCs did not exceed the baseline (0.55+/-0.04 0.56+/-0.14, 0.58+/-0.38) with CIs above the null hypothesis or p-values below 0.05 (Supplementary Figure 2 A-J, Supplementary Table 1). For the tumor types PAAD, OV and PRAD, the AUROCs ranged from 0.58+/-0.22 to 0.6+/-0.09 to 0.76+/-0.22. Together, these data demonstrate that DL can predict HRD status from histology images alone in several tumor types.

• Download figure
• Open in new tab

Figure 2: Comparison of Area under the receiving operating curve (AUROC) for internal and tumor wise external validation experiment models.

Boxplot displaying the distribution for the AUROC for (A) internal 5-fold cross-validation experiment of The Cancer Genome Atlas (TCGA) and tumor-wise external validation on the Clinical Proteomic Tumor Analysis Consortium (CPTAC); (B) AUROCs for the cross-cancer external validation experiment of the TCGA breast cancer cohort (TCGA-BRCA) on the TCGA and CPTAC cohort. The horizontal line indicates the median, whereas each box represents the interquartile range (IQR) between the first and third quartiles. The whiskers extend from the box to the minimum and maximum values, considering 1.5 times the IQR. Abbreviations: BRCA=breast cancer; CRC=colorectal cancer; GBM=glioblastoma; LIHC=liver cancer; LUAD=lung adenocarcinoma; LUSC/LSCC=lung squamous cell carcinoma; OV=ovarian cancer; PAAD/PDA=pancreatic adenocarcinoma; PRAD=prostate adenocarcinoma; UCEC=endometrial cancer

HRD is predictable from H&E histology with attMIL in an independent test set

A step that is germane to the successful development of DL models is external validation with WSIs from patient cohorts which are completely independent from the training set (45). Hence, for our external validation experiments, we deployed the classification models obtained from the cross-validation training on TCGA to analyze cohorts from the CPTAC dataset corresponding to the same cancer type. External validation cohorts in CPTAC were available for endometrial cancer (UCEC), pancreatic cancer (PDA), lung adenocarcinoma (LUAD), and lung squamous cell carcinoma (LSCC). In these external validation experiments, we noted that the prediction performance was highercompared to internal validation experiments. Once again, the best performance was obtained in UCEC, with an AUROC of 0.93+/-0.07, p=0.01. In LUAD, the performance increased in the external validation, yielding an AUROC of 0.73+/-0.11 and a significant p-value of 0.03. In the case of PAAD/PDA, where the internal validation was unsuccessful (internal validation AUROC 0.58+/-0.22), the external validation resulted in an improved AUROC reaching 0.81+/-0.14, albeit with a p-value of 0.07. Meanwhile, in LUSC/LSCC, no improvement in performance was observed in the external validation set compared to the internal training set (AUROC 0.57+/-0.01, p=0.23, Figure 2A, Supplementary Figure 2 K-N). Together, these data show that DL-based classifiers of HRD status generalize beyond the training cohort.

A HRD classifier trained on BRCA detects HRD across various types of cancer

Lastly, we aimed to investigate if HRD-related morphological features in one cancer type can help to predict HRD status in another cancer type. This would point to a shared set of morphological features across cancer types, potentially allowing us to develop a pan-cancer pathology-based prediction system for HRD status. To test this, we applied our trained HRD classifiers in a cross-cancer experimental design. We used the breast cancer cohort TCGA-BRCA to train the HRD classification model because this cohort had the highest number of patients. Subsequently, we deployed this model on all other cohorts obtained from the TCGA and CPTAC datasets. Surprisingly, the BRCA-based model was able to significantly predict HRD from non-BRCA tissue in UCEC, PRAD and PAAD. For those three cohorts, the external deployment of a BRCA-based model resulted in higher prediction AUROCs than the respective internal validation experiments, reaching AUROCs of 0.70+/-0.02, p<0.001 in TCGA-UCEC, 0.84+/-0.07 and p=0.004 in TCGA-PRAD 0.67+/-0.03, p=0.2 in TCGA-PAAD, 0.87+/-0.1 p=0.05 in CPTAC-UCEC and 0.65+/-0.02 p=0.26 in CPTAC-PDA, respectively (Figure 2B). In the tumor types LUAD and OV, the AUROCs remained with 0.62+/-0.03 for TCGA-LUAD, 0.66+/-0.06 for CPTAC-LUAD and 0.61+/-0.03 in TCGA-OV in a similar range to the internal validation results (Supplementary Figure 3A-M). Together, these data show that a classifier trained on breast cancer can predict HRD status from histology in other tumor types, indicating a shared “HRD morphology” between tumor types.

Molecular and histomorphological characterization of TCGA-BRCA HRD-H and HRD-L cases

Finally, we investigated which molecular and morphological patterns were associated with ground truth and DL-predicted HRD status. We used the TCGA-BRCA cohort to analyze this in detail, as this was the largest cohort. We observed that in the HRD-H subgroup, 45% were classified as basal-like breast cancers, 11% as HER2-enriched, 15% as Luminal A, and 26% as Luminal B. In contrast, only 7% of the cases in the HRD-L subgroup were basal-like, 7% were HER2-enriched, 64% were Luminal A, and 18% were Luminal B (Figure 3A) (46). Within our predicted groups, we observed a similar distribution among the BRCA subtypes (Figure 3B).

• Download figure
• Open in new tab

Figure 3: Molecular Characterization of The Cancer Genome Atlas breast cancer (TCGA-BRCA) cohort.

(A) Distribution of breast cancer subtypes for the Homologous Recombination deficiency high (HRD-H) and low (HRD-L) ground truth subgroups. (B) Distribution of the breast cancer subtypes for the HRD-H and HRD-L Deep Learning (DL) predicted subgroups. (C) Alteration Frequency for several genes of the HRD-H and HRD-L ground truth subgroups. (D) Alteration Frequency for several genes of the HRD-H and HRD-L within cohort internal results prediction subgroups. (E) Grouped Boxplots comparing the Homologous Recombination Deficiency high (HRD-H) prediction scores with the mutational status (mutated=MUT, wildtype=WT) for the somatic and germline alterations of the BRCA1/2 genes. The central line represents the median value, while the box ranges between the first and third quartile (IQR) and the whiskers extend to the lowest and highest values within 1.5 times the IQR. The y-axis represents the Deep Learning (DL) HRD-H prediction values. An independent t–test was performed to calculate the p-values: ns: p <= 1.00e+00 *: 1.00e-02 < p <= 5.00e-02 **: 1.00e-03 < p <= 1.00e-02 ***: 1.00e-04 < p <= 1.00e-03

To reassure that our model predicts HRD detached from phenotypic differences of estrogen receptor negative (ER-) vs. ER-positive (ER+) breast cancers we calculated the receiving operating curve (ROC) and precision recall curve (PRC) for the subgroups: ER+/HER2+, ER+/HER2-, ER-/HER2+, ER-/HER2-, indicating that HRD was also predictable with AUROCs of 0.66+/-0.3, 0.8+/-0.09, 0.72+/-0.43 and 0.62+/-0.11 (Supplementary Figure 4A-H). Our analysis of the mutational landscape of both HRD-H and HRD-L ground truth revealed that TP53 had the highest alteration frequency with 67% in the HRD-H ground truth group, significantly higher than 20% in the HRD-L group, following alterations in the TTN (26% vs. 14%) gene. In contrast, the most enriched alterations in the HRD-L group were observed in the genes PIK3CA (39%) followed by CDH1 (16%), GATA3 (14%) and MAP3K1 (11%), whereas the prevalences in the HRD-H group of PIK3CA, CDH1, GATA3 and MAP3K1 were 19%, 2%, 6% and 1%, respectively (Figure 3C). For the HRD-H prediction subgroup alteration frequencies for TP53, were significantly higher at 77% (Figure 3D). Such divergences were not as noticeable in the HRD-L prediction group. These findings suggest that there are notable differences in alteration frequencies between the two subgroups, which are consistent across both the ground truth and prediction data. Moreover, we compared the HRD-H prediction score to the alteration status of somatic and germline mutations in the BRCA1/2 genes, whereupon we saw that there was a significant difference between the mutant and wild-type cases for BRCA1 germline and BRCA2 somatic mutations (Figure 3E). Methylation data indicated that the HRD-H group had most of its methylation alterations in the N-shore portion of the BRCA1 promoter region, whereas those in the HRD-L group were mainly located in the S-shore portion (Supplementary Figure 4I). Lastly, we proceeded to investigate the histomorphological patterns associated with the presence of HRD through whole slide prediction heatmaps in CPTAC-UCEC (Figure 4A-C). Our findings revealed that high grade, fibrosis, hemorrhage and lymphocytic infiltration are consistent features predictive of HRD across various tumor types, as shown in Figure 4 for BRCA and UCEC, particularly in the top predicted HRD-H tiles for the top three patients. Fibrosis was observed in HRD-positive cases, particularly in BRCA (Figure 4D). Moreover, hemorrhagic necrosis especially adjacent to tumor tissue and tumor stroma was consistently observed as highly predictive areas in the true HRD-H cases across various cancer types. (Supplementary Figure 5, 6). In summary, these data show that known HRD morphology characteristics were found in our DL based top predicted HRD-H patients.

• Download figure
• Open in new tab

Figure 4: Visualization of predicted Homologous Recombination Deficiency high (HRD-H) tumor samples.

(A) Whole slide image (WSI) of an HRD-H predicted patient (ID: C3L-00358-21) from the Clinical Proteomic Tumor Analysis Consortium (CPTAC) endometrial cancer (UCEC) cohort with magnification. (B) Attention heatmap for the same patient with magnification. (C) Classification Heatmap for the same patient with magnification. (D) Top predicted tiles for top three homologous recombination deficiency high (HRD-H) patients in The Cancer Genome Atlas (TCGA) breast cancer (BRCA). (E) Top predicted tiles for three HRD-H patients in the CPTAC-UCEC cohort.