Tristan Lazard

Guillaume Bataillon

Peter Naylor

Tatiana Popova

François-Clément Bidard

Dominique Stoppa-Lyonnet

Marc-Henri Stern

Etienne Decencière

Thomas Walter
email: thomas.walter@mines-paristech.fr

Anne Vincent Salomon
email: anne.salomon@curie.fr

Deep Learning identifies new morphological patterns of Homologous Recombination Deficiency in luminal breast cancers from whole slide images

Keywords: homologous recombination deficiency (HRD), breast cancer, computational pathology, artificial intelligence, deep learning, prediction, visualization, bias, molecular subtype, whole slide images, interpretability

Homologous Recombination DNA-repair deficiency (HRD) is a well-recognized marker of platinum-salt and PARP inhibitor chemotherapies in ovarian and breast cancers (BC). Causing high genomic instability, HRD is currently determined by BRCA1/2 sequencing or by genomic signatures, but its morphological manifestation is not well understood. Deep Learning (DL) is a powerful machine learning technique that has been recently shown to be capable of predicting genomic signatures from stained tissue slides. However, DL is known to be sensitive to dataset biases and lacks interpretability. Here, we present and evaluate a strategy to control for biases in retrospective cohorts. We train a deep-learning model to predict the HRD in a controlled cohort with unprecedented accuracy (AUC: 0.86) and we develop a new visualization technique that allows for automatic extraction of new morphological features related to HRD. We analyze in detail the extracted morphological patterns that open new hypotheses on the phenotypic impact of HRD.

The advent of Deep Learning has revolutionized biomedical image analysis and in particular digital pathology. Traditionally, the majority of methods developed in this field were dedicated to computer-aided diagnosis, where the objective is to partially automatize human interpretation of slides, in order to help pathologists in their diagnosis task, e.g. the detection of mitoses [START_REF] Veta | Assessment of algorithms for mitosis detection in breast cancer histopathology images[END_REF], or the identification of metastatic axillary lymph nodes [START_REF] Bejnordi | Diagnostic Assessment of Deep Learning Algorithms for Detection of Lymph Node Metastases in Women With Breast Cancer[END_REF][START_REF] Campanella | Clinical-grade computational pathology using weakly supervised deep learning on whole slide images[END_REF]. Beyond the automatization of manual inspection, Deep Learning has also been successfully applied to predict patient variables, such as outcome [START_REF] Mobadersany | Predicting cancer outcomes from histology and genomics using convolutional networks[END_REF], and to predict molecular features, such as gene mutations [START_REF] Kather | Pan-cancer image-based detection of clinically actionable genetic alterations[END_REF][START_REF] Coudray | Classification and mutation prediction from non-small cell lung cancer histopathology images using deep learning[END_REF], expression levels [START_REF] Schmauch | A deep learning model to predict RNA-Seq expression of tumours from whole slide images[END_REF] or genetic signatures [START_REF] Kather | Pan-cancer image-based detection of clinically actionable genetic alterations[END_REF][START_REF] Diao | Human-interpretable image features derived from densely mapped cancer pathology slides predict diverse molecular phenotypes[END_REF]. Despite these results of unprecedented quality, one of the major drawbacks of Deep Learning algorithms is their black-box character: because the features are automatically extracted, it is difficult to know how a decision was made. This has two major consequences: first, it is difficult to identify potential confounders, i.e. variables that correlate with the output due to the composition of the data set and that are predicted instead of the intended output variable. Second, even in the absence of statistical artifacts, understanding how the decision was generated in the first place can point to interesting mechanistic hypotheses and to patterns in the image that have so far been overlooked.

One way to overcome the latter problem is to use hand-crafted biologically meaningful features [START_REF] Diao | Human-interpretable image features derived from densely mapped cancer pathology slides predict diverse molecular phenotypes[END_REF]. This however requires an extraordinary effort in terms of annotation. Here, we take a conceptually different approach. Instead of working in a pan-cancer setting on a large number of signatures, we concentrate on one single medically highly relevant signature in one cancer type on a controlled data set, where we can correct for potential biases. In order to understand how the Deep Learning decision is generated and which morphological patterns are related to the output variable, we propose a novel visualization and interpretation technique that paves the way to "machine teaching", i.e. a data driven approach to identify phenotypic patterns related to genomic signatures, thus pointing to new mechanistic hypotheses.

In order to demonstrate the power of this strategy, we focus on predicting Homologous Recombination Deficiency (HRD) in Breast Cancer (BC). Worldwide, 2.1 million women are newly diagnosed per year with BC which is a leading cause of cancer-related death. Improvement of metastatic breast cancer treatment is therefore of highest priority. BC is a heterogeneous disease with four major molecular classes (luminal A, B, HER2 enriched, and triple-negative breast cancer [TNBC]) benefiting from different therapeutic approaches.

If early BC patients have an overall survival of 70 to 80%, metastatic disease is incurable with a short duration of survival [START_REF] Deluche | Contemporary outcomes of metastatic breast cancer among 22,000 women from the multicentre ESME cohort 2008-2016[END_REF]. Homologous Recombination (HR) is a major and high-fidelity repair pathway of DNA double-strand breaks. Its deficiency, HRD, results in high genomic instability [START_REF] Miller | ESMO recommendations on predictive biomarker testing for homologous recombination deficiency and PARP inhibitor benefit in ovarian cancer[END_REF] and occurs through diverse mechanisms, including germline or acquired mutations in DNA repair genes, most frequently BRCA1, BRCA2 or PALB2, or through epigenetic alterations of BRCA1 or RAD51C. Importantly, HRD leads to a high sensitivity to polyADP-ribose polymerase inhibitors (PARPi) in vitro [START_REF] Bryant | Specific killing of BRCA2-deficient tumours with inhibitors of poly(ADP-ribose) polymerase[END_REF][START_REF] Farmer | Targeting the DNA repair defect in BRCA mutant cells as a therapeutic strategy[END_REF]. PARPi have been shown to improve metastatic breast cancer progression free survival [START_REF] Tung | TBCRC 048: Phase II Study of Olaparib for Metastatic Breast Cancer and Mutations in Homologous Recombination-Related Genes[END_REF][START_REF] Tutt | Adjuvant Olaparib for Patients with BRCA1-or BRCA2-Mutated Breast Cancer[END_REF]. Several methods have been developed to detect HRD, including genomic instability profiling, mutational signatures, or integrating structural and mutational signatures [START_REF] Popova | Ploidy and Large-Scale Genomic Instability Consistently Identify Basal-like Breast Carcinomas with BRCA1/2 Inactivation[END_REF][START_REF] Birkbak | Telomeric Allelic Imbalance Indicates Defective DNA Repair and Sensitivity to DNA-Damaging Agents[END_REF][START_REF] Abkevich | Patterns of genomic loss of heterozygosity predict homologous recombination repair defects in epithelial ovarian cancer[END_REF][START_REF] Polak | A mutational signature reveals alterations underlying deficient homologous recombination repair in breast cancer[END_REF][START_REF] Davies | HRDetect is a predictor of BRCA1 and BRCA2 deficiency based on mutational signatures[END_REF]. However, HRD is currently diagnosed in clinical practice by DNA repair genes sequencing and genomic instability patterns (genomic scar) such as the LST signature [START_REF] Popova | Ploidy and Large-Scale Genomic Instability Consistently Identify Basal-like Breast Carcinomas with BRCA1/2 Inactivation[END_REF] or the HRD MyChoice® CDx test (Myriad Genetics). BRCA1 and BRCA2 mutations are known predictive markers for response to PARPi [START_REF] Miller | ESMO recommendations on predictive biomarker testing for homologous recombination deficiency and PARP inhibitor benefit in ovarian cancer[END_REF][START_REF] Tutt | Adjuvant Olaparib for Patients with BRCA1-or BRCA2-Mutated Breast Cancer[END_REF] and platinum salt [START_REF] Tutt | Carboplatin in BRCA1/2-mutated and triple-negative breast cancer BRCAness subgroups: the TNT Trial[END_REF] and the somatic HRD has been more recently recognized as a predictive marker for PARPi in ovarian [START_REF] Miller | ESMO recommendations on predictive biomarker testing for homologous recombination deficiency and PARP inhibitor benefit in ovarian cancer[END_REF] and breast cancer [START_REF] Chopra | Homologous recombination DNA repair deficiency and PARP inhibition activity in primary triple negative breast cancer[END_REF]. But neither a specific routinely assessed phenotype nor a morphological pattern indicates the presence of HRD. The majority of hereditary BRCA1 cancers are TNBC and up to 60-69% of sporadic TNBC harbor a genomic profile of HRD [START_REF] Popova | Ploidy and Large-Scale Genomic Instability Consistently Identify Basal-like Breast Carcinomas with BRCA1/2 Inactivation[END_REF][START_REF] Alexandrov | Signatures of mutational processes in human cancer[END_REF][START_REF] Chopra | Homologous recombination DNA repair deficiency and PARP inhibition activity in primary triple negative breast cancer[END_REF] . However, the majority of hereditary BRCA2 cancers are luminal [START_REF] Lakhani | The pathology of familial breast cancer: predictive value of immunohistochemical markers estrogen receptor, progesterone receptor, HER-2, and p53 in patients with mutations in BRCA1 and BRCA2[END_REF] and HRD also exists in sporadic luminal B [START_REF] Chopra | Homologous recombination DNA repair deficiency and PARP inhibition activity in primary triple negative breast cancer[END_REF][START_REF] Manié | Genomic hallmarks of homologous recombination deficiency in invasive breast carcinomas[END_REF] or in HER2 tumors [START_REF] Ferrari | A whole-genome sequence and transcriptome perspective on HER2-positive breast cancers[END_REF][START_REF] Turner | Signatures of DNA-Repair Deficiencies in Breast Cancer[END_REF]. Of note, germline or sporadic alterations of BRCA harbor indistinguishable genomic alterations in triple-negative or in luminal tumors [START_REF] Manié | Genomic hallmarks of homologous recombination deficiency in invasive breast carcinomas: Genomic hallmarks of homologous recombination defect[END_REF][START_REF] Holstege | BRCA1-mutated and basal-like breast cancers have similar aCGH profiles and a high incidence of protein truncating TP53 mutations[END_REF]. In that context, a reliable and accurate test is mandatory to select patients for PARPi and platinum salt treatments. Whereas the screening for germline or somatic BRCA1 and BRCA2 mutations is feasible for all TNBC (18% of all BC cases), it represents a real challenge in clinical practice if extended to all luminal B tumors (35% of all BC). The recent results of the Olympia trial emphasize the need for an efficient method of screening for BRCA1 and 2 mutations across all breast cancer phenotypes [START_REF] Tutt | Adjuvant Olaparib for Patients with BRCA1-or BRCA2-Mutated Breast Cancer[END_REF]. This strategy moreover does not identify the whole diversity of genetic causes of HRD.

In this study, we present an image-based approach to predict HR status from Whole Slide Images (WSI) stained with Hematoxylin Eosin (HE) using deep learning, from a large retrospective series of luminal and triple-negative breast carcinomas with a genomically defined HR status, from a single cancer center. In particular, we show that careful correction for potential biases is essential for such studies and demonstrate the relevance of a correction strategy. Finally, we develop a novel interpretation algorithm that allows the visualization of decisive previously unknown patterns and leads to new hypotheses on disease-relevant genotype-phenotype relationships.

Results

A Deep Learning architecture to predict HRD from Whole Slide Images (WSI)

We scanned the most representative HE stained tissue section of the surgical resections specimens of breast cancer from 715 patients with known HR status. The series was composed of 309 Homologous Recombination Proficient (HRP) tumors and 406 Homologous Recombination Deficient tumors (Supplementary Table 1).

Due to their enormous size, analysis of WSI typically relies on the Multiple Instance Learning (MIL) paradigm [START_REF] Ilse | Attention-based Deep Multiple Instance Learning[END_REF][START_REF] Amores | Multiple instance classification: Review, taxonomy and comparative study[END_REF][START_REF] Maron | A Framework for Multiple-Instance Learning[END_REF][START_REF] Courtiol | Classification and disease localization in histopathology using only global labels: a weakly supervised approach[END_REF]. MIL techniques only require slide-level annotations and share the overall architecture (see Figure 1) consisting of 4 main steps: tiling and encoding, tile-scoring, aggregation, and decision.

[ Figure 1: From WSI to prediction. Four major components are used in this end-to-end pipeline. First, the WSI ( ) are tiled, the tissue parts are automatically selected, and the resulting tiles are embedded into a low-dimensional space (block 1). The embedded tiles are then scored through the attention module ( 2). An aggregation module outputs the slide level vector representative (3) that is finally fed to a decision module (4) that outputs the final prediction. When training, the binary cross-entropy loss between the ground truth and the prediction is computed and back-propagated to update the parameters of the modules. Both the decision module and the attention module are multi-layer perceptrons, the encoder is a ResNet18 and the aggregation module consists of a weighted sum of the tiles, the weights being the attention scores. ]

The WSI was divided into tile images (dimension: 224x224 pixels) arranged in a grid.

Background tiles are removed, tissue tiles are encoded into a feature vector. Instead of using representations trained on natural image databases and unlike most studies in this domain, we used the self-supervised technique Momentum Contrast (MoCo [START_REF] He | Momentum Contrast for Unsupervised Visual Representation Learning[END_REF]; see Methods). This method consists in training a Neural Network to recognize images after transformations, such as geometric transformations, noise addition and color changes. By choosing the kind and strength of transformations, we can impose invariance classes, i.e. variations in the input that do not result in different representations. The feature vector of each tile was then mapped to a score by a neural network. The slide representation was obtained by the sum of the individual tile representations, weighted by the learned attention scores [START_REF] Ilse | Attention-based Deep Multiple Instance Learning[END_REF]. Finally, the slide representation was classified by the decision module (see Figure 1). We optimized hyperparameters by a systematic random search strategy (see Methods). For hyperparameter setting and performance estimation, we used nested 5-fold cross-validation, which allows us to obtain realistic performance estimations. All reported performance results are averaged over 5 independent test folds (see Methods).

HRD prediction with correction for potential biases

We applied this method to predict HRD from the WSI in the TCGA cohort, and obtained results (AUC=0.71; see Figure 2.c) in line with previous reports [START_REF] Diao | Human-interpretable image features derived from densely mapped cancer pathology slides predict diverse molecular phenotypes[END_REF][START_REF] Valieris | Deep Learning Predicts Underlying Features on Pathology Images with Therapeutic Relevance for Breast and Gastric Cancer[END_REF][START_REF] Kather | Deep learning can predict microsatellite instability directly from histology in gastrointestinal cancer[END_REF]. While the TCGA is an invaluable resource for pan-cancer studies in genomics and histopathology, it is often seen rather as a starting point whose results need to be corroborated by other cohorts [START_REF] Kleppe | Designing deep learning studies in cancer diagnostics[END_REF]. We then devised a sampling strategy that mitigates biasing during training. Bias mitigation is an increasingly important line of research in machine learning. For instance, it is a well-known problem in training predictive models for functional magnetic resonance imaging (fMRI) data, where the age of the patient has been shown to be an important confounder [START_REF] Varoquaux | Assessing and tuning brain decoders: cross-validation, caveats, and guidelines[END_REF]. While several techniques for bias mitigation exist [START_REF] Zhao | Training confounder-free deep learning models for medical applications[END_REF][START_REF] Zhao | Men Also Like Shopping: Reducing Gender Bias Amplification using Corpus-level Constraints[END_REF][START_REF] Adeli | Representation Learning with Statistical Independence to Mitigate Bias[END_REF][START_REF] Wang | Balanced Datasets Are Not Enough: Estimating and Mitigating Gender Bias in Deep Image Representations[END_REF], a recent comparison [START_REF] Wang | Towards Fairness in Visual Recognition: Effective Strategies for Bias Mitigation[END_REF] indicates that strategic sampling is the method of choice if the distribution is not too In addition to these technical confounders, we identified the molecular subtype of the tumor to be a potential biological confounder. Successful correction of this biological confounder in the TCGA (see figure 2.b) led however to a dramatic drop in performance (AUC=0.63). This result suggests that NN trained on the TCGA for HRD prediction without stratification or bias correction, might actually predict to some extent the molecular subtype, which is also a predictable variable (AUC=0.89). This shows that the molecular subtype is indeed a biological confounder. In our in-house data set, we decided to build a subtype-specific NN that specifically predicts HRD for luminal BC, instead of applying bias mitigation. The reason for this decision was three-fold: first, we argued that a dataset focusing on only one Attention scores are often used for visualization in the field of pathology [START_REF] Campanella | Clinical-grade computational pathology using weakly supervised deep learning on whole slide images[END_REF][START_REF] Lu | Data-efficient and weakly supervised computational pathology on whole-slide images[END_REF][START_REF] Dehaene | Self-Supervision Closes the Gap Between Weak and Strong Supervision in Histology[END_REF][START_REF] Courtiol | Deep learning-based classification of mesothelioma improves prediction of patient outcome[END_REF], either in the form of heatmaps in order to localize the origin of the relevant signals or in the form of galleries of tiles of interest (tiles with highest attention scores).

However, attention scores do not per se extract the tiles that are related to a certain output variable; they just reflect that the tile is to be taken into consideration in the decision. In particular in the case of genetic signatures, where we would expect that the output variables can be related to several morphological patterns, analyzing only the attention scores might be limited. In the case of HRD prediction, this intuition is corroborated by the results presented in Figure 3: selecting the tiles with the highest attention scores does not allow to identify tiles related to HRD or HRP, respectively. Moreover, a low-dimensional representation, obtained by Uniform Manifold Approximation and Projection (UMAP) (47) of the tile descriptors does not show any grouping that seems to relate to the output variable.

[Figure 3 Our approach thus suggests that these tissue phenotypes are hallmarks of HRD. In contrast to triple-negative tumors that have been described as high-grade, rich in TILs, with pushing margins [START_REF] Rakha | Patho-biological aspects of basal-like breast cancer[END_REF], no specific pathological patterns of luminal breast gBRCA 1 or 2 cancers were identified except a high grade with a frequent absence of tubules formation and pushing margins. However, these features did not allow a robust identification of HRD luminal tumors in clinical practice [START_REF] Stratton | Pathology of familial breast cancer: differences between breast cancers in carriers of BRCA1 or BRCA2 mutations and sporadic cases[END_REF][START_REF] Lakhani | The Pathology of Familial Breast Cancer: Histological Features of Cancers in Families Not Attributable to Mutations in BRCA1 or BRCA2[END_REF][START_REF] Bane | BRCA2 Mutation-associated Breast Cancers Exhibit a Distinguishing Phenotype Based on Morphology and Molecular Profiles From Tissue Microarrays[END_REF]. Here, thanks to our visualization tool and our unbiased deep-learning analysis, we identified new features linked to HRD in luminal breast cancers.

In order to test this analysis result, the TIL density and nuclear grade were evaluated for each tumor of the in-house dataset by an expert pathologist. As predicted by our algorithm,

Discussion

In this study, we set out to predict the HR status in breast cancer from H&E stained WSI and to analyze the phenotypic patterns related to HRD. The prediction of HRD is an important challenge in clinical practice. The use of PARP inhibitors for breast cancer patients was initiated for metastatic TNBC patients with germline mutations of BRCA1 or BRCA2.

However, BRCA2, as well as PALB2 and a minority of BRCA1 cancer patients, develop luminal tumors. The necessity of predicting HRD is therefore not limited to TNBC, but extends also to luminal BC. On the other hand, luminal BCs represent a far more frequent group than TNBC. For this reason, systematic screening of HR gene alterations for luminal cancers will be problematic and in many countries even infeasible due to both economic and logistic issues. Therefore, preselection of patients with a high probability of being HR deficient by analysis of WSI is a cost-efficient strategy that has so far only been hampered by the lack of knowledge about HRD specific morphological patterns in luminals. Indeed, only high grades and to a lower extent pushing margins have previously been reported to be associated with HRD. In this context, the identification of HRD from WSI by deep learning and the identification of related morphological patterns could both facilitate the preselection of breast cancers for molecular determination of HRD, which is particularly important for luminal cancers.

The TCGA provides a precious data set to train models for the prediction of genetic signatures from H&E data [START_REF] Kather | Pan-cancer image-based detection of clinically actionable genetic alterations[END_REF][START_REF] Diao | Human-interpretable image features derived from densely mapped cancer pathology slides predict diverse molecular phenotypes[END_REF]. While we obtained promising results for the prediction of HRD on the TCGA dataset in line with previous reports, we found that this result was partly due to the fact that the molecular subtype acts as a biological confounder. This was particularly problematic as we wanted to investigate the morphological signature of HRD. Of note, the existence of biological and technical confounders is presumably not limited to HRD prediction, but may concern many genetic signatures. The use of carefully curated data sets where technical and biological confounders can be controlled for, is thus an important step in investigating the predictability of genetic signatures, as well as the identification of their morphological counterparts.

In most cases, such in-house datasets also contain technical and biological biases, due to the long period during which the dataset is acquired. This motivated us to propose a method to mitigate bias in Computational Pathology workflows, based on strategic sampling. Such strategies are already used in other fields of medical imaging but have so far -to the best of our knowledge -not been used in Computational Pathology. We have shown that this approach can successfully mitigate or even eliminate bias. However, the method is limited by the number of variables we can correct for, as well as by the class imbalance it can handle.

In some cases, stratification might therefore be preferable. In any case, it is important to be aware of the confounding variables in the data set, whose presence can lead to false expectations and misinterpretation. For this reason, we expect proper treatment of such variables to become a standard in the field.

While bias correction on the TCGA led to a drop in AUC to 0.63, we found that HRD was predictable in our in-house data set of 715 BC patients with an AUC of 0.83. While homogeneous datasets do not reflect the variability between centers and thus limit direct applicability of the trained networks, they allow for controlled feasibility studies, which now need to be complemented by multicentric studies. In addition, we will validate this algorithm in a prospective neoadjuvant clinical trial for which patients' HRD status will be assessed with MyChoice® CDx test (Myriad).

Homogeneous datasets are well suited for the identification of phenotypic patterns linked to disease patterns, even in cases where no such patterns are known a priori, such as in the case for HRD. In order to identify a phenotypic signature related to an output variable (here HRD), we can either use biologically meaningful encodings, also known as human interpretable features (HIF), and infer the most relevant features by analyzing the weights in the predictive model ( 8), or we can turn to network introspection. The HIF approach relies on detailed and exhaustive annotations of a large number of WSI. For instance, (8) leverage annotations provided by hundreds of pathologists consisting of hundreds of thousands of manual cell and tissue classifications. Here, we provide a new network introspection scheme, relying on the powerful MoCo encodings, trained without supervision directly on histopathology data, and a decision-based tile selection, that allows us to automatically cluster tiles and to relate these clusters to the output variable.

Interestingly, while our approach confirms the recently published finding that necrosis is a hallmark of HRD ( 8) and identifies morphological features common to HRD in TNBC and luminal BC, such as necrosis, high density in TILs and high nuclear anisokaryosis [START_REF] Rakha | Patho-biological aspects of basal-like breast cancer[END_REF], it also points to more specific patterns that have so far been overlooked. For instance, we found tiles enriched in carcinomatous cells with clear cytoplasm suggesting activation of specific metabolic processes in these cells. Second, we find intra-tumoral laminated fibrosis as an HRD related pattern. This suggests the hypothesis that cancer-associated fibroblast Similar to what we have shown here with respect to HRD, the visualization framework we have developed is versatile and can in principle be applied in the context of other genetic signatures. Because the algorithm is fully automated, using the MIL algorithm and its visualization method can constitute a useful tool for the discovery of morphological features related to the predicted genetic signatures. This has the potential to generate new biological hypotheses on the phenotypic impact of these genetic disorders. In order to maximize the benefit for the scientific community, we release the code to train MIL models on WSIs and to create morphological maps as well as tile trajectories publicly and free of charge, and provide detailed documentation.

Altogether, this study provides new and versatile tools for the prediction and phenotypic dissection of genetic signatures from histopathology data. Application to luminal breast cancers allowed us to shed light on the phenotypic consequences of homologous recombination deficiency, and might provide a tool with the potential to impact breast cancer patient care.

Methods Data

In-house dataset (Institut Curie). We retrospectively retrieved a series of 715 patients with HE slides of surgical resections specimens of untreated breast cancer and a genomically known HR status (Supplementary Table 1 All patients have been treated and followed at the Institut Curie between 1995 and 2020.

The patient agreed for the use of tumor samples from their surgical resection specimens for research according to the law. Ethical approval from the Institutional Review Board (Institut Curie breast cancer study group N°DATA190031) was obtained for the use of all specimens.

Clinical data have been retrieved from the Institut Curie electronic medical records and saved using Research electronic data capture (REDCap) tools hosted at the Institut Curie.

Public dataset (TCGA)

. This public dataset is composed of 815 WSI of breast cancer fixed in formalin (FFPE) and stained in H&E. They are available at https://portal.gdc.cancer.gov/ .

Low-resolution WSI, WSI containing artifacts such as pen marks, tissue-folds and blurred WSI were removed. The final data-set encompasses 691 WSIs. The HR status of the corresponding tumors was obtained using the LST genomic signature [START_REF] Popova | Ploidy and Large-Scale Genomic Instability Consistently Identify Basal-like Breast Carcinomas with BRCA1/2 Inactivation[END_REF].

Architecture and optimization parameters. Training and performance evaluation are done in a 5-fold nested cross-validation framework.

Each dataset is split into 5 independent folds. For each of these folds, a validation set is randomly sampled in the complementary 4/5th. A model is trained on the remaining dataset ( = ⅘ * ⅘ th of the total dataset). This process is repeated 10 times for each test fold, then the best model is selected according to its validation performances, and finally tested on its test set.

Each test and validation set preserves the stratification of the whole dataset with respect to the target variable as well as the confounding variables in case we correct for them. The final performance estimation of the model is the performance averaged over the 5 test performances. During inference time, all the tiles of each WSI are processed Strategic sampling is used both for balancing the training dataset with respect to the output variable ( in the binary case) and to correct for biases ( ).

If is a given WSI sampled from the dataset, then and are respectively the target value and the bias value of . We note the total number of slides in the dataset labeled with , same for . is the total number of slides with label value and bias value .

For achieving both, we sample the WSIs in each batch in a distribution under which and for all i.

That is : for each , .

Strategic sampling is performed on the fly when building the batches.

When correction for several confounder simultaneously, and

, a new confounder variable is created that takes values in all the pairs combinations of and :

. We then apply strategic sampling to correct for C.

Bias score. The bias score of a confounder variable is the average mutual information between and the predicted class ( ) estimated in a dataset in which the mutual information between and the target is zero.

We sample 30 sub-dataset in the test set using strategic sampling such that and average over them to get the bias score :

The mutual information between two finite random variables B and Q taking value respectively in and measures the dependency between B and Q and is defined by : .

  Figure 2.c). As the cohort was generated over 25 years, two experimental variables ( , ,

  imbalanced. Strategic sampling aims at ensuring that irrespectively of the composition of the training set, each batch presented to the neural network is composed of roughly the same number of samples for each value combination of output and confounding variable. Correcting for and resulted in a 4-fold reduction of the bias-score in comparison to the uncorrected model and a slightly lower accuracy (AUC=0.86, Figure 2.c).

[Figure 2 :

 2 molecular subtype was more likely to reveal the underlying patterns exclusively related to HRD; second, HRD prediction in luminal BC is of particular importance for clinical practice, as very few morphological patterns are known to be related to HRD in luminal BC, the most frequent breast cancer phenotype, and third, the relatively low number of TNBC in our data set made strategic sampling on three confounding variables challenging. Therefore, we composed a dataset containing only luminal BC and setting both technical confounders, leading us to keep 259 BC WSI (191 HRD tumors, 66 HRP tumors). We obtained a good, albeit slightly lower performance of this bias-corrected NN (AUC=0.83, Figure2.c; Figure 2.d). The trained model carefully freed from both technical and biological biases was then used for the identification of morphological patterns described in the next section. Bias corrections and prediction performances. a-b, estimation of the bias score of two technical confounder ( , ) and one biological confounder ( ) for the Curie Dataset (a) and the bias score of the confounder for the TCGA dataset (b) for different correction strategies. A Mann-Whitney-Wilcoxon test two-sided with Bonferroni correction is performed for each pair of correction strategies. ns: non-significant p>0.05, *: p<0.05, **: p<0.01, ***: p<1e-3, ****: p<1e-4 c-d, performance results. Name of each model indicates the origin of its training set. Indices indicate the correction applied through strategic sampling ( has been debiased with respect to ). corresponds to the model trained on a subset containing only luminal tumors. c: ROC curve of the models trained on the Curie dataset correcting for technical bias ( ) or for technical biases and ( ). d: summary tables of performance metrics. AUC : Area Under The (receiver operating characteristics) Curve, : balanced accuracy.]Visualization reveals HRD specific tissue patternsIn order to understand which phenotypic patterns are related to HRD on the WSI, we turned to visualization techniques for NN. The used MIL framework is equipped with an inherent visualization mechanism: the second module of the algorithm, the tile-scoring module, is in fact an attention module that assigns to each tile an attention score that determines how much a given tile will contribute to the slide representation (and thus to the decision).

  : I. Attention-based visualization does not discriminate between HRD and HRP. I. a, Mechanism of the attention-based visualization. The attention score of a tile is used as a direct proxy of its importance in the prediction of the WSI. I.b, UMAP projection of the highest attention ranked tiles of the Curie WSIs classified as HRP (orange crosses) and HRD (blue circles). I.c Randomly sampled tiles among the HRP and HRD tiles. The tiles are located in the tumor, however, neither clear clusters nor visual differences are present between HRD and HRP tiles. II. Decision-based visualization. II.a, Mechanism of the decision-based visualization. 1, Each tile in the whole dataset is scored by the attention module. 2, The best scoring tiles are selected as candidate tiles. 3, The selected tiles are presented to the decision module, the probability of each of these tiles being HRD or HRP (yellow or green) is kept. 4, Finally, the K tiles with maximal probability for either HRD / HRP are selected. II.b, Morphological map of the HR status. Each dot is the UMAP projection of a tile extracted by the decision-based visualization method. Crosses are tiles with high probability of being HRD; circles are tiles maximizing HRP prediction. Each cluster has been linked to a morphological phenotype by two expert pathologists. We identified six different morphological phenotypes associated with the HRD and two associated with the HRP. The exhibited tiles have been randomly sampled among each cluster. II.c Pathological interpretation of the clusters presented in II.b.Given these limitations, we propose a new visualization protocol that allows us to extract the tiles that are directly associated with a particular slide-level label. As the slide representation is the weighted sum of the tile representations, we applied the decision module, specifically trained to classify slide representations between HRD and HRP, to the individual tile representations. This gives us a score for each tile that can be interpreted as the (tile) probability of being HRD or HRP (see Methods for details). Selecting the tiles with the highest posterior probability for HRD and HRP respectively, and projecting the tile representations of this selection to a low dimensional space leads to the emergence of distinct clusters corresponding to different tumor tissue patterns with a clear relation to HRD or HRP and therefore providing a morphological map of HRD (Figure4). Two expert pathologists labeled these clusters. The HRD signal relied on several clusters: HRD tumors present a high tumor cell density, with a high nucleus/cytoplasm ratio and conspicuous nucleoli. They also show regions of hemorrhagic suffusion associated with necrotic tissue. In the stroma, the HRD signal revealed the presence of striking laminated fibrosis and as expected relied on high Tumor-Infiltrating Lymphocytes (TILs) content. Lastly, one large cluster contained a continuum of several phenotypes, namely adipose tissue intermingled with scattered and clear tumor cells, histiocytes, and plasma cells. In contrast, the HRP signal was mostly carried by one cluster characterized by low tumor cell density, the cells being moderately atypical and tumor cell nests separated from the stroma by clear spaces. Notably, it included a few invasive lobular carcinomas (all the tiles per cluster are available in the supplementary figures 1-11).

[Figure 4 :

 4 TILs and nuclear grade were positively associated with the HR status of the tumor in the luminal subset (mean TILs HRD: 29, mean TILs HRP: 17, t-test-pvalue : 0.017; mean nuclear grade HRD: 2.7, mean nuclear grade HRP: 2.3, Xi2-pvalue : 1.2 e-6 ).Our NN works with different internal representations. While the tile representations provided by MoCo permit the emergence of phenotypic similarity clusters (Figure4), internal representations closer to the decision module encode information relevant for HRD. The representation in the penultimate layer can therefore be interpreted as encoding "HRD-ness" of the tiles. Figure4illustrates a low dimensional representation of this HRD-ness for the same tiles as those present in Figure4, where point color represents the HRD-score (tile probability to be classified as HRD). From there, we have extracted two tile trajectories, going from low HRD-ness to high HRD-ness. The magenta trajectory illustrates the successive visual changes corresponding to an increase in tumor cells or inflammatory cells density (from low-density tiles to high-density tiles with large nuclei, nuclear atypia and infiltrative lymphocytes). The blue trajectory shows conversely a decrease in tumor cells density replaced successively by an inflammatory reaction and apoptotic cells, loose fibrosis and hemorrhagic suffusion associated with necrosis. These different trajectories illustrate the manifestations of HRD and show the pleiotropic character of the induced phenotypes. Moreover, the highlighted gradation of these phenotypes opens the path to a possible reading grid of WSIs for pathologists. Illustration of 2 Phenotypic HRD-ness trajectories. a: UMAP projection of the HR status specific representation of the meaningful tiles relative to the HRD. HRD-ness is the score given to each tile by the HRD output neuron. Two tile trajectories have been extracted (blue and magenta) starting from the same low HRD-ness region, each leading to a different high HRD-ness region. b, c: Tiles sampled along each of the trajectories. They are ordered from low HRD-ness to high HRD-ness and read from left to right and from up to bottom. b. Magenta trajectory, toward densely cellular tumors or inflammatory cells. c: Blue trajectory, toward fibroinflammatory tumor changes and hemorrhagic suffusions.

(

  CAF) within the stroma of HRD luminal tumors may play a role in the viability and fate of tumor cells. Furthermore, the presence of adipose tissue within the tumor suggests first a different tumor cell density and second a specific balance between CAF and adipocytes in the context of a luminal HRD tumor. The molecular mechanisms achieving these patterns remain to be determined by in vitro models.

  Hyperparameters have been set thanks to a random search evaluated through 5-fold nested cross-validation. The benchmark task is the prediction of the molecular class of the TCGA WSIs. Both the decision module and the tile-scoring module are multilayer perceptrons with batch normalization[START_REF] Ioffe | Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift[END_REF] after each hidden layer. The decision module has 3 hidden layers of 512 neurons, the tile-scoring module has 1 hidden layer of 256 neurons.Dropout has been fixed at 0.4, the optimizer is ADAM[START_REF] Kingma | ADAM: A Method for Stochastic Optimization[END_REF] with a learning rate of 3e-3. A batch consists of 16 samples of WSI. A sample of WSI corresponds to a uniform sampling of 300 of its composing tiles. In fact, we observed that this uniform subsampling of the WSIs regularized training as well as diminishes its computational workload. Finally, training is performed during 200 epochs.

Acknowledgements

The authors thank Hélène Guenon, Saida Sahiri Martial Caly and Laure Annette for their help in retrieving the H&E slides and their technical expertise. GB was supported by a Fondation Curie grant. TL was supported by a Q-Life PhD fellowship (Q-life ANR-17-CONV-0005). Furthermore, this work was supported by the French government under management of Agence Nationale de la Recherche as part of the ''Investissements d'avenir'' program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute).

Code availability: The code that allows to reproduce this study and perform new ones is available at ****************************** Supplementaries: Supplementaries are available at https://drive.google.com/file/d/1e1gEU8KfKl6LMq3lODWV3CJor9BgT_pE/view?usp=sharing

Because this average is non-negative even when is not biased (but tends to zero when is large), we compute the bias score of an unbiased model predicting the target variable with the same accuracy as the tested model. It serves as a control.

Learning MoCo representations. For learning MoCo-v2 representation we used the MoCo repository available at https://github.com/facebookresearch/moco. We randomly used the following transformations: Gaussian blur, crop and resize, color jitter, grayscale, horizontal and vertical symmetries, and finally a color augmentation in the Hematoxylin and Eosin specific space [START_REF] Ruifrok | Quantification of histochemical staining by color deconvolution[END_REF].

The training dataset is composed of 5.3e6 images of size 224x 224 pixels, or half the Curie dataset at magnification 10x.

We used a Resnet18 and trained it for 60 epochs on 4 GPU Nvidia Tesla V100 SXM2 32 Go. We used the SGD optimizer with a momentum of 0.9, a weight decay of 1e-4 and a learning rate of 3e-3. We used a cosine scheduler with warm restart on the learning rate.

Visualization methods.

The model used to extract the visualizations has been trained on the luminal subset of the Curie dataset (259 WSI). To benefit from the biggest dataset possible, the model has been trained on the whole dataset, without using early stopping nor testing, during 200 epochs.

To generate the attention-based visualisation, the highest ranked tile with respect to the attention score is extracted, for each WSI. The selected tiles are then labeled according to the label of their WSI of origin.

Concerning the decision-based visualisation, for each WSI the 300 highest ranked tiles with respect to the attention score are selected. Amond this pool of tiles, the 2000 highest ranking tiles with respect to the posterior probability for HRD and HRP are selected. In order to promote diversity in the extracted images, no more than 20 tiles per slide can be selected.

Computation resources. All computations have been done on the GENCI HPC cluster of

Jean-Zay.

Authors contributions

AVS and GB initiated the project. AVS, GB, FCB and DSL generated the patient cohort. GB reviewed all the slides. TP and MHS performed the genomic analyses. TL, ED and TW designed the AI and statistical methods. TL and PN developed the software. TL performed the analysis, under the supervision of TW and ED. AVS and GB interpreted the morphological patterns. AVS, GB, TL, ED and TW discussed methods, results and design choices. TL prepared the figures. TL, TW and AVS wrote the manuscript.


Mélanie Lubrano Di Scandalea

Tristan Lazard

Guillaume Balezo

Yaëlle Bellahsen-Harrar

Cécile Badoual

Sylvain Berlemont

Thomas Walter

Automatic grading of cervical biopsies by combining full and self-supervision

In computational pathology, the application of Deep Learning to the analysis of Whole Slide Images (WSI) has provided results of unprecedented quality. Due to their enormous size, WSIs have to be split into small images (tiles) which are first encoded and whose representations are then agglomerated in order to solve prediction tasks, such as prognosis or treatment response. The choice of the encoding strategy plays a key role in such algorithms. Current approaches include the use of encodings trained on unrelated data sources, full supervision or self-supervision. In particular, self-supervised learning (SSL) offers a great opportunity to exploit all the unlabelled data available. However, it often requires large computational resources and can be challenging to train. On the other end of the spectrum, fully-supervised methods make use of valuable prior knowledge about the data but involve a costly amount of expert time. This paper proposes a framework to reconcile SSL and full supervision and measures the tradeoff between long SSL training and annotation effort, showing that a combination of both has the potential to substantially increase performance. On a recently organized challenge on grading Cervical Biopsies, we show that our mixed supervision scheme reaches high performance (weighted accuracy (WA): 0.945), outperforming both SSL (WA: 0.927) and transfer learning from ImageNet (WA: 0.877). We further provide insights and guidelines to train a clinically impactful classifier with a limited expert and/or computational workload budget. We expect that the combination of full and self-supervision is an interesting strategy for many tasks in computational pathology and will be widely adopted by the field.

Introduction

Recent advances in slide digitization have led to increased interest in Artificial Intelligence (AI) applications for histopathology. The development of AI models could help reduce pathologists' workloads, limit subjectivity and help contributing to medical discoveries. Deep learning models can now match pathologist performance for many tasks: diagnostic, detection of mitoses [START_REF] Veta | Assessment of algorithms for mitosis detection in breast cancer histopathology images[END_REF], prediction of gene mutations [START_REF] Coudray | Classification and mutation prediction from non-small cell lung cancer histopathology images using deep learning[END_REF][START_REF] Nikolas Kather | Pan-cancer image-based detection of clinically actionable genetic alterations[END_REF] or genetic signatures [START_REF] Nikolas Kather | Pan-cancer image-based detection of clinically actionable genetic alterations[END_REF][START_REF] Diao | Human-interpretable image features derived from densely mapped cancer pathology slides predict diverse molecular phenotypes[END_REF][START_REF] Lazard | Deep Learning identifies new morphological patterns of Homologous Recombination Deficiency in luminal breast cancers from whole slide images[END_REF], cancer subtyping [START_REF] Coudray | Classification and mutation prediction from non-small cell lung cancer histopathology images using deep learning[END_REF] and more. One of the applications, automated diagnosis from Whole Slide Images (WSIs), induces two main challenges: first, WSIs are very high-resolution and, because of memory constraint, cannot be fed directly into traditional neural networks. Second, expert annotations are laborious to attain, costly and prone to subjectivity. The most popular methods today rely on Multiple Instance Learning (MIL), which frames the problem as a bag classification task. WSIs are split into small workable images (tiles), which are processed separately. Features from each of the individual tiles are extracted and then aggregated to classify the WSI. The extraction of these tiles' specific representation is crucial to the downstream WSI classification task. One common approach consists of initializing the feature extractor with pre-trained weights on ImageNet, a natural image dataset. This technique allows one to extract generic features that are powerful, but that do not lie within the histopathological domain. Different strategies have been developed to extract these tile encodings taking advantage of the available data and their respective level of supervision. A first strategy aims to learn tile features with full supervision [START_REF] Babak | and and the CAMELYON16 Consortium. Diagnostic Assessment of Deep Learning Algorithms for Detection of Lymph Node Metastases in Women With Breast Cancer[END_REF]. To create a supervised dataset, one or several experts manually review tiles and sort them into meaningful classes (preferably related to the downstream task of classifying the WSIs). Even though experts' annotations can bring powerful prior knowledge to the model, this technique often requires large quantities of annotations. A second strategy consists of learning tile representations through self-supervision. It leverages the unannotated data by training a convolutional neural network on a pretext task. It has proven its efficacy [START_REF] Saillard | Identification of pancreatic adenocarcinoma molecular subtypes on histology slides using deep learning models[END_REF][START_REF] Lu | Data-efficient and weakly supervised computational pathology on whole-slide images[END_REF] and even its superiority to the fully supervised scheme [START_REF] Dehaene | Self-Supervision Closes the Gap Between Weak and Strong Supervision in Histology[END_REF]. However, this approach has a non-negligible computational cost, as training necessitates around 1000 hours of computation on a standard GPU [START_REF] Dehaene | Self-Supervision Closes the Gap Between Weak and Strong Supervision in Histology[END_REF]. Moreover, it is not guaranteed that the obtained encodings are most relevant for the prediction task we are trying to solve. Techniques from both sides of the supervision spectrum have proven to bring important benefits for relevant feature extraction. Combining them could allow us to benefit from the best of both worlds. In this work, in addition to proposing a joint-optimization process mixing self, full and weak supervision (Figure 1), we measure the trade-off in performance between the number of annotations and the computational cost of training a self-supervised model. We thus provide guidelines to train a clinically impactful classifier with a limited budget in expert and/or computational workload.

2

Attention weights transfer

Unlabelled Tiles

Labelled Tiles

Labelled Slides

Step 3: Whole Slide Classification

Step 2: Fine-tuning

Related Work

Mixed Supervision Medical data is often limited. For this reason, one might want to take advantage of all the available data even if annotations might not be homogeneous and even though they might be difficult to exploit because multiple levels of supervision are available. For instance, whole slide images are often associated with one global label (weak supervision), they can contain millions of unlabelled tiles (no supervision), but, as a pathologist reviews the slides and performs a diagnostic, it is almost effortless for them to mark the region of interest that signs the corresponding 3 diagnostic (strong supervision). AI applications have usually been dichotomized between supervised and unsupervised methods, spoiling the potential of combining several types of annotations. For this reason, mixing supervision for medical images analysis has gained interest in past years [START_REF] Huang | Rectifying Supporting Regions With Mixed and Active Supervision for Rib Fracture Recognition[END_REF][START_REF] Li | Thoracic Disease Identification and Localization with Limited Supervision[END_REF], 2021]. For instance, in [START_REF] Mlynarski | Deep learning with mixed supervision for brain tumor segmentation[END_REF] the author showed that combining global labels and local annotations by training in a multi-task setting, the capacities of the model to segment brain tumors on Magnetic Resonance Images were improved.

In [START_REF] Tourniaire | Attention-based Multiple Instance Learning with Mixed Supervision on the Camelyon16 Dataset[END_REF], the author introduced a mixed supervision framework for metastasis detection building on the CLAM [START_REF] Lu | Data-efficient and weakly supervised computational pathology on whole-slide images[END_REF] architecture. CLAM is a variant of the popular attention based MIL [START_REF] Ilse | Attention-based Deep Multiple Instance Learning[END_REF] with 2 extensions: first, in order to make the method applicable in a multi-class setting, class-specific attention scores are learned and applied. Second, the last layer of the tile encoding network is trained to also predict the top and bottom attention scores, thus mimicking tile-level annotations. In [START_REF] Tourniaire | Attention-based Multiple Instance Learning with Mixed Supervision on the Camelyon16 Dataset[END_REF], the authors highlight the limitations of this instance-classification approach and propose to leverage a low number of fully annotated slides to train the attention mechanism. In a second step, they propose to turn to a standard MIL training (using only slide-level annotations). Even with few annotated slides, this approach allows to boost classification performance. However, there are also some limitations. First, the method relies on exhaustive annotation of selected slides: for the annotated slides, all the key regions are annotated pixel-wise. Second, due to the CLAM architecture, the approach only fine-tunes a single dense layer downstream the pre-trained feature extractor. Third, the algorithm has been designed for an application case in which the slide and tile labels coincide (tumour presence). This however is not always the case: when predicting genetic signatures, grades or treatment responses, it is unclear how tile and slide level annotations relate to each other. In this article, we propose to overcome these limitations. We propose to combine self-supervised learning with supervision prior to training the MIL network. We thus start from more powerful encodings, that are not only capable of solving the pretext task of self-supervised learning, but also the medical classification task that comes with the annotated tiles. Consequently, this method does not require full-slide annotations, optimizes the full tile encoding network and does not come with any constraint regarding the relationship between tile and slide level annotations.

3 Materials and Method

Dataset and Problem Setting

The Tissue Net Challenge [DrivenData] organized in 2020, the Société Française de Pathologie (SFP) and the Health Data Hub aimed at developing methods to automatically grade lesions of the uterine cervix in four classes according to their severity. The training dataset for the challenge was made up of biopsy samples from female uterine cervix, focusing on squamous lesions (Figure 2). These lesions are often benign but can also be qualified as low grade or high grade depending on the 4 risk of invasion of the underlying conjunctive tissue and evolution into carcinomas. The grade of the lesions depends on the proportion of squamous epithelium affected by dysplastic criteria. Lowgrade squamous intraepithelial lesions (LSIL) are defined as having a dysplastic criteria involving less than one third of the thickness of the epithelium. High-grade squamous intraepithelial lesions (HSIL) indicate a greater proportion of the epithelium composed of undifferentiated basal cells with abnormalities. Carcinoma is diagnosed when abnormal epithelial cells invade the underlying conjunctive tissue. The class of a WSI was determined by the highest lesion's grade present on it.

Fully Supervised Dataset

5926 annotated Regions of Interest (ROIs) of fixed size 300x300 micrometers were provided. Each ROI had roughly the same size as a tile at 10x magnification and were labeled by the severity of the lesion it contained: "Normal" (0) if tissue was normal, (1) LSIL or (2) HSIL if it presented precancerous lesions that could have malignant potential and (3) invasive squamous carcinoma (Table 1).

Weakly Supervised Dataset

The dataset was composed of 1015 WSIs acquired from 20 different centers in France at an average resolution of 0.234 +/-0.0086 mpp (40X). The slide resolution varied slightly due to the multicentric provenance of the data. The class of the WSI corresponded to the class of the most severe lesions it contained (grade from 0 to 3 also). All the native WSI formats were converted to pyramidal TIFF (Tagged Image File Format). Both the WSI-level and tile-level labels have been attributed by a consortium of expert pathologists (Table 1).

Misclassification Costs

Misclassification errors do not lead to equally serious consequences (i.e predicting a benign lesion if it is cancerous is more serious than predicting a LSIL instead of a HSIL). Accordingly, a panel of pathologists established a grading of each of these errors i.e they attributed to each pair of possible outcome (i, j) ∈ {0, 1, 2, 3} 2 a severity score 0 C i,j 1 (Table 2)

The metric used in the challenge to evaluate and rank the submissions is computed from the average of these misclassification costs.

More precisely, if we name P (S) the prediction of a slide S labelled l(S), the challenge metric M W A is:

M W A = 1 N S (1 -C l(S),P (S) ) (1)
with N the number of samples.

The problem is thus framed as a cost-sensitive classification problem, and, to our knowledge, all the winning solutions took awareness of this cost in their training procedure. In Multiple Instance Learning, we are given sets of samples

B k = {x i |i = 1 . . . N k }, also called bags.
The annotation y k we are given refers only to the bags and not the individual samples. We assume however, that such tile-level labels exist in principle, but that we just do not have access to them.

6

The strategy is to first map each tile x i to its encoding z i , which is then mapped to a scalar value a i , often referred to as attention score. The tile representations z i and attention scores a i are then agglomerated to build the slide representation s k which is then further processed by a neural network. The agglomeration can be based on tile selection [START_REF] Campanella | Clinical-grade computational pathology using weakly supervised deep learning on whole slide images[END_REF], Courtiol et al., 2020], or on an attention mechanism [START_REF] Ilse | Attention-based Deep Multiple Instance Learning[END_REF], which is today the most widely used strategy.

Self-Supervised Learning

Self-supervised learning provides a framework to train neural networks without human supervision. The main goal of self-supervised learning is to learn to extract efficient features with inputs and labels derived from the data itself using a pretext task. Many self-supervised approaches are based on contrastive learning in the feature space. SimCLR, a simple framework relying on data augmentation was introduced in [START_REF] Chen | A Simple Framework for Contrastive Learning of Visual Representations[END_REF]. Powerful feature representations are learned by maximizing agreement between differently augmented views of the same data point via a contrastive loss applied in the feature space.

An image is transformed through random data augmentations into two new images. They are then embedded using the feature extractor. The two features vectors (z i and z j ) are mapped with a projection head (dense layers) to obtain final vectors h i and h j . The feature extractor and projection head are trained to maximize agreement using the contrastive loss. Positive pairs consist of the two augmented views of the same image, the other 2(n -1) views play the role of negative samples. The loss function (NT-Xent) for a positive pair (i, j) is defined as:

L SSL = -log exp(sim(h i , h j )/τ 2n k=1 1 k =i exp(sim(h i , h j ) τ (2)
Where sim(u, v) = u T v ||u||•||v|| , the cosine similarity, 1 k =i∈(0,1) determines if k = i and τ is a parameter. After convergence, the projection head is discarded and the pretrained feature extractor can be used for subsequent tasks.

Cost-Sensitive Training

Instead of the traditional cross-entropy loss we used a cost-aware classification loss, the Smooth-One-Sided Regression Loss L SOSR . First introduced to train SVMs in [START_REF] Tu | One-sided Support Vector Regression for Multiclass Costsensitive Classification[END_REF], this objective function was smoothed and adapted for backpropagation in deep networks in [START_REF] Chung | Cost-Aware Pre-Training for Multiclass Cost-Sensitive Deep Learning[END_REF]. When using this loss, the network is trained to predict the class-specific risk rather than a posterior probability; the decision function chooses the class minimizing this risk. The SOSR loss is defined as follows:

L SOSR = i j ln(1 + exp(2 i,j • (ĉ i -C i,j ))) (3) 7
With 2 i,j = -1 i =j + 1 i=j , ĉi the i-th coordinate of the network output and C the error table.

Mixed Supervision

To be tractable, training of attention-MIL architectures requires freezing the feature extractor weights. While SSL allows the feature extractor to build meaningful representations [START_REF] Saillard | Identification of pancreatic adenocarcinoma molecular subtypes on histology slides using deep learning models[END_REF][START_REF] Dehaene | Self-Supervision Closes the Gap Between Weak and Strong Supervision in Histology[END_REF], they are not specialized to the actual classification problems we try to solve. Several studies have shown that such SSL models benefit from fine-tuning specific to the downstream task [START_REF] Chen | A Simple Framework for Contrastive Learning of Visual Representations[END_REF] We therefore added a training step to leverage the tile-level annotation and fine-tune the selfsupervised model. However, as the final WSI classification task is not identical to the tile classification task, we suspect that fine-tuning solely on the tile classification task may over-specialize the feature extractor and thus sacrifice the generalizability of SSL (and for this reason ultimately also degrading the WSI classification performances). To avoid this, we developed a training process that optimizes the self-supervised and tile-classification objectives jointly. Two different heads, plugged before the final classification layer, are used to compute both loss functions L SSL and L SOSR The final objective L is then:

L = βL SSL + (1 -β)L SOSR (4)
where β is a hyperparameter that has to be tuned. Here, we found β = 0.3 (see Supplementary).

Understanding the Feature Extractor with Activation Maximization

To further understand the features learned by the different pre-training policies (ImageNet, supervised, SSL and mixed), we used Activation Maximization (AM) to visualize extracted features and provide an explicit illustration of the specificity learned. Methods to generate pseudo-images maximizing a feature activation have been introduced in [START_REF] Dumitru Erhan | Visualizing Higher-Layer Features of a Deep Network[END_REF]. This technique consists in synthesizing the images that will maximize one feature activation. It is summarized as follow [START_REF] Nguyen | Understanding Neural Networks via Feature Visualization: A survey[END_REF]]: If we consider a trained classifier with set of parameters θ that map an input image x ∈ R h×w×c , (h and w are the height and width and c the number of channels) to a probability distribution over the classes, we can formulate the following optimization problem:

x * = arg max x (σ l i (θ, x)) (5)
where σ l i (θ, x) is the activation of the neuron i in a given layer l of the classifier. This formulation being a non-convex problem, local maximum can be found by gradient ascent, using the following 8 update step:

x t+1 = x t + ∂σ l i (θ, x) ∂x t (6)
The optimization process starts with a randomly initialized image. After a few steps, it generates an image which can help to understand what information is being captured by the feature. As we try to visualize meaningful representations of the features, some regularization steps are applied to the random noise input (random crop and rotations to generate more stable visualization, details can be found in Supplementary Materials). To generate filter visualization within the HE space, we transformed the RGB random image to HE input thanks to color deconvolution [START_REF] Ruifrok | Quantification of histochemical staining by color deconvolution[END_REF]. This preprocessing allowed to generate images with histology-like colors when converted back to the RGB space. To select the most meaningful features for each class, we trained a Lasso classifier without bias to classify the extracted feature vectors into the four classes of the dataset for the four pre-training policies. The feature vectors for each tile were first normalized and divided element-wise by the vector of features' standard deviation across all the tiles. The L1 regularization factor λ was set to 0.01. Details about Lasso training can be found in Supplementary Materials. Contribution scores for each feature were therefore derived from the weights of the Lasso linear classifier: negative weights were removed and remaining positive weights were divided by their sum to obtain contribution scores [0,1]. By filtering out the negative weights, the contribution score corresponds to the proportion of attribution among the features positively correlated to a class, and allows to select feature capturing semantic information related to the class, leaving out those containing information for other classes.

6 Experimental Setting

WSI Preprocessing

Preprocessing on a downsampled version of the WSIs was applied to select only tissue area and non-overlapping tiles of 224x224 pixels were extracted at a resolution of 1 mpp. (Details in Supplementary Materials)

Data Splits for Cross-Validation

To measure the performances of our models we performed 3-fold cross-validation for all our training settings. Because the annotated tiles used in our joint-optimization step were directly extracted from the slides themselves, we carefully split the tiles such that tiles in different folds were guaranteed to originate from different slides. The split divided the slides and tiles into a training set, a validation set and a test set. All subsequent performance results are then reported as the average and standard deviation of the performance results on each of these 3 test folds. 9

Feature Extractor Pre-Training

The feature extractor is initialized with pre-trained weights obtained with three distinct supervision policies: fully supervised, self-supervised or a mix of supervision. These three policies rely on the fine-tuning of a DenseNet121 [START_REF] Huang | Densely Connected Convolutional Networks[END_REF], pretrained on ImageNet. The fully-supervised archicture is fine-tuned solely on the tile classification task. The SSL architecture is derived from SimCLR framework and is trained on an unlabeled dataset of 1 million tiles extracted from the slides. Finally for the mixed-supervised architecture, a supervised branch is added to the previous SSL network and trained using the mixed objective function (see Fig. 1 and Eq. 4) on the fully supervised dataset. Technical details of these three training settings are available in the supplementary material.

Whole Slide Classification

After tiling the slides, the frozen feature extractor (DenseNet121) was applied to extract meaningful representations from the tiles. This feature extractor was initialized sequentially with the pre-trained weights mentioned above and generated as many sets of features. These bags of features were then used to train the Attention-MIL model with SOSR loss applied slide-wise. (Supplementary Materials).

Feature Visualization

To select the most relevant features, we trained an unbiased linear model on the feature vectors extracted from the annotated tiles. The feature vectors were standardised. The weights of the linear model were used to determine which features were the most impactful for each class. Feature visualizations were generated for the selected features and for each set of pre-trained weights. We extracted the tiles expressing the most of these features by selecting the feature vectors with the higher activation for the concerned feature. Implementation details are provided in Supplementary Materials.

Results

Self-Supervised Fine-Tuning

We saved the checkpoints of the self-supervised feature extraction model at each epoch of training, allowing us to investigate the amount of time needed to reach good WSI classification performances. We computed the embeddings of the whole dataset with each of the checkpoints and trained a WSI classifier from them. Figure 3 reports the performances of WSI classification models for each of these checkpoints. SSL training led to a higher Weighted Accuracy than using ImageNet weights after 3 epochs and resulted in a gain of +4.8% after 100 epochs. Interestingly, as little as 6 epochs 10 of training are enough to gain 4% of Weighted Accuracy: a significant boost in performance is possible with 50 GPU-hours of training. We then observe a small increase in performance until the 100th epochs.

Pre-Training Policy Comparison

To compare the weights obtained with the various supervision levels, we ran a 3-fold cross-validation on the WS classification task and summarized the results in table 3. The results indicate that the SSL pre-training substantially improves the WSI classification performance. In contrast, we see that initializing the feature extractor with fully-supervised weights gives an equivalent or poorer performance than any other initialization. SSL pre-training allows us to extract rich features that are generic, yet still relevant to the dataset (unlike ImageNet). On the other hand, fully supervised features are probably too specific and seem to not represent the full diversity of the image data. The joint-optimization process manages to balance out generic and specialized features without neutralizing them: mixing the supervision levels brings significant improvements (+2%) to the performance, leading to a Weighted Accuracy of 0.945. We additionally compared the benefits introduced by the cost-sensitive loss (Eq. 3) with the crossentropy loss. Our results show that with ImageNet weights the SOSR loss improves the Weighted Accuracy by 1% and the accuracy by 3%. In conclusion, the combination of the SSL pre-trained model, its fully supervised fine-tuning, and the cost-sensitive loss leads to a notable improvement of 8 Weighted Accuracy points over the baseline MIL-imagenet model. (Chen, Kornblith, Swersky, et al. 2020) that showed that an SSL model is up to 10x more label efficient than a supervised one. However, for the 100-epoch checkpoint, we observe an improvement of 2 points of the Weighted Accuracy when using 100% of the tile annotations. Moreover, fine-tuning the models by mixed supervision with too few annotations (15%) leads to a slight drop in WSI classification performances. Finally, we see a diminution of the standard deviations across splits for the different pre-training policies, showing better stability for longer SSL training and more annotations.

We draw different conclusions from these observations:

• In this context, it is always better to pre-train the feature extractor with SSL rather than only invest in annotations.

• The supervised fine-tuning needs enough annotations to bring an improvement to the WSI classification task. We can note however that even when considering the 100% annotation settings, the supervised dataset (approx. 5000 images) is still rather small in comparison to traditional image datasets.

• A full SSL training is mandatory to leverage this small amount of supervised data.

Features Visualisations

We generated the pseudo-images of the most important features for each class and each pre-training policy and extracted the related tiles. The Figure 5 displays the most important features along with 12 Table 4: Relationship between self-supervision and full-supervision -Study on the performance improvement on WS classification for different proportion of labelled data versus different training time of SSL the tiles activating each feature the most for the class "Normal" (0). Although interpretation of such pseudo images must be treated carefully, we notice that the features obtained with SSL, supervised and mixed training are indubitably more specialized to histological data than those obtained with ImageNet. Some histological patterns, such as nuclei, squamous cells or basal layers are clearly identifiable in the generated images. The extracted tiles are strongly correlated with class-specific biomarkers. Feature e represents a normal squamous maturation, i.e. a layer of uniform and rounded basal cells, with slightly larger and bluer nuclei than mature cells. We can also observe several layers of mature cells with small nuclei and moderately abundant cytoplasm (pink halo around), equidistant from each other. Features c and d highlight clouds of small regular and rounded nuclei (benign cytological signs). Feature g and h are characteristic of squamous cells (polygonal shapes, stratified organization lying on a straight basal layer). Interestingly, features extracted with the supervised method (g, h) manage to sketch a normal epithelium with high resemblance, the features are more precise. On the other hand, features extracted with SSL (c, d) highlight true benign criteria but do not entirely summarize a normal epithelium (no basal maturation). The mixed model displays both, suggesting that mixed supervision highlights pathologically relevant patterns to a larger extent than the other regimes [WHO, 2020]. We can also note by looking at the real tiles that while features from ImageNet (a, b), SSL (c, d) and the supervised model (g, h) focus on the upper half of the cervix epithelium, it appears that features from the mixed supervision model (e, f ) are focusing on the lower half which is known to be the relevant region for discrimination between class Normal (0) and Low Grade (1) (abnormal cells are constricted to the lower third of the epithelium).

In Figure 4 we can further identify class-related biomarkers for dysplasia and carcinoma grade. Tiles with visible koilocytes (cells with a white halo around the nucleus) have been extracted from the top features for Low Grade class. Koilocytes are symptomatic of infection by Human Papillomavirus and are a key element for this diagnosis (almost always responsible for precancerous lesions in the cervix, [WHO, 2020]). High Grade (2) generated image represents disorganised cells with a high nuclear-to-cytoplasmic ratio, marked variations in size and shape and loss of polarity. For the class "Carcinoma" (3), we observe irregular clusters of cohesive cells with very atypical nuclei, separated by a fibrous texture that can be identified as stroma reaction. All these criteria have been identified in [WHO, 2020] as key elements for diagnosis of dysplasia and invasive carcinoma.

13

In Figure 6, we observe that features extracted from ImageNet and SSL models are diverse, in particular, features extracted from SSL reflect rich tissue phenotypes which correlates to their generic capacities of image representations. On the other hand, features extracted with supervised and mixed methods are more redundant. We additionally observe in Figure 6

Discussion

In pathology , expert annotations are usually hard to obtain. However, we are often in a situation where a small proportion of labeled annotation exists but not in sufficient quantities to support fully supervised techniques. Yet, even in small quantities, expert annotations carry meaningful information that one could use to enforce biological context to deep learning models and make sure that networks learn appropriate patterns. On the other hand, self-supervised methods have proven their efficacy to extract generic features in the histopathological domain and their usefulness for 16 downstream supervision tasks, even in the absence of massive ground truth data. Methods capable of reconciling self-supervision with strong supervision can therefore be useful and open the door to better performances.

In this paper, we presented a way to inject the fine-grained tile level information by fine-tuning the feature extractor with a joint optimization process. This process allowed to mix self-supervised learning features with tile classification ones and helped the downstream WSI classification task.

We our method to the TissueNet Challenge, challenge for the automatic grading of cervix cancer, that provided annotations at the slide and tile level, thus representing an appropriate use case to validate our method of mixed supervision. We also propose in this study insights and guidelines for the training of a WSI classifier in the presence of tile annotations.

First, we showed that SSL is always beneficial to our downstream WSI classification tasks. Fine-tuning pre-trained weights with SSL for only 50 hours brings a 4% improvement over WSI classification weighted accuracy, and near to 5% when fine-tuning for longer (100 epochs).

Second, a small set of annotated tiles can bring benefit to the WSI classification task, up to 2% of weighted accuracy for a supervised dataset of around 5000 images. Such a set of tiles can be obtained easily by asking the pathologist to select a few ROIs that guided his decision while labeling the WSIs, which can be achieved without a strong time commitment. However this boost in performance can be reached only if the feature extractor is pre-trained with SSL, and for sufficiently long: SSL unlocks the supervised fine-tuning benefits.

To further understand the differences between the range of supervision used to extract tile features, we conducted qualitative analysis on features visualizations by activation maximization and observed that features obtained from SSL, supervised or mixed trainings were more relevant for histological tasks and that class-discriminative patterns were indeed identified by the model. We also observed that supervised training on the tiles alone led to much less diversity in the features extracted by the model than the ones obtained with SSL.

The scope of this study contains by design three limitations. First, SSL models were trained by fine-tuning already pre-trained weights on imagenet. This may explain the rapid convergence and boost in performance observed; however it may also underestimate this boost if the SSL models were trained from scratch. We did not compare SSL trained from scratch and fine-tuned SSL, and left it to future work.

Second, all the conclusions reached are conditioned by the fact that we do not fine-tune the feature extractor network during the WSI classification training. Keeping these weights frozen, and even pre-computing the tiles representations brings a large computational benefit (both in memory and speed of computations), but prevents the feature extractor from specializing during the WSI classification training. Third, the tendency observed in table 4 of better performances correlated with larger numbers of annotations is modest and would require more annotations to validate it. Finally, our method can be improved in several ways. First, SimCLR, was a pioneer method in self-supervised learning architecture and has proven to be efficient but it suffers from high performance drop when decreasing the batch size [START_REF] Chen | A Simple Framework for Contrastive Learning of Visual Representations[END_REF]. Other SSL models have been developed to alleviate this limitation. MoCo [START_REF] He | Momentum Contrast for Unsupervised Visual Representation Learning[END_REF] actually propose a momentum mechanism allowing optimal performances even without large batch size and therefore, numerous available parallel GPUs. Other models like VICReg [START_REF] Bardes | VICReg: Variance-Invariance-Covariance Regularization for Self-Supervised Learning[END_REF] proposed techniques to maximize the variance between the features and therefore limit their redundancies. It will be interesting to benchmark these SSL variants with respect to their impact on WSI classification accuracy and feature interpretability.

To conclude, we present a method that provides an interesting alternative to using full supervision, pre-training on unrelated data sets or self-supervision. We convincingly show that the learned feature representations are both leading to higher performance and providing more intermediate features that are more adapted to the problem and point to relevant cell and tissue phenotypes. We expect that the mixed supervision will be adopted by the field and lead to better models.

Figure 1 :

 1 Figure 1: Mixed Supervision Process: a)A self-supervised model (SimCLR) is trained on unlabelled tiles extracted from the slides. Feature extractor and contrastive layer weights are transferred to the joint-optimization architecture b) Joint-optimization model is trained on the labeled tiles of the dataset. The feature extractor weights are transferred to the WS classification model. c) WS classification model is trained on the 1015 whole slide images.

Figure 2 :

 2 Figure 2: Illustration of Uterine cervix dysplasia -[National Cancer Institute, 2011]

Figure 3 :

 3 Figure 3: Weighted Accuracy evolution -Weighted Accuracy evolution on WS classification task with respect to the number of epochs of SSL training

  that feature visualisation from the mixed model picture realistic histopathological patterns specific to the class. Visualisation for other classes are available in Supplementary Materials.

Figure 4 :Figure 5 :

 45 Figure 4: Feature comparison per class -The top row displays the top filter for the Mixed Supervised model for each class. The bottom row displays the tile expressing the feature the most

Figure 6 :

 6 Figure 6: Feature Diversity for the class "Carcinoma" (3) (top 5 features) -Class "Normal" (0) and top 10 features in Supplementary Materials

Table 1 :

 1 Dataset Summary

	Classes	Number of Slides Number of Tiles
	0 (Normal)	270	1923
	1 (Low Grade)	288	1405
	2 (High Grade)	238	1368
	3 (Carcinoma)	219	1230
	Total	1015	5926

Table 2 :

 2 Weighted Accuracy Error Table -Error table to ponderate misclassification according to their gap with the ground truth

	4 Proposed Architecture

4.1 Multiple Instance Learning and Attention

Table 3 :

 3 Pre-training policies -Performances summary7.3 Number of Annotations vs Number of EpochsWe have seen that both SSL and supervised pre-training bring together an improvement in the WSI classification task. To further investigate the relationship between these two supervision regimes, we trained models with only some of the fully supervised annotations (15, 65, 100%) on top of intermediate SSL checkpoints. Results are reported in table 4. It appears that without SSL pre-training (or with too few epochs of training), the supervised finetuning does not bring additional improvement in WSI classification. This is in line with the work of Chen et al.

	11

Acknowledgments

The authors thank Etienne Decencière for the thoughful discussions that helped the project. ML was supported by a CIFRE PhD fellowship founded by KEEN EYE, Paris, France and ANRT (CIFRE 2019(CIFRE /1905)). TL was supported by a Q-Life PhD fellowship (Q-life ANR-17-CONV-0005). Furthermore, this work was supported by the French government under management of Agence Nationale de la Recherche as part of the "Investissements d'avenir" program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute).