Arthur Imbert

Wei Ouyang

Adham Safieddine

Emeline Coleno

Christophe Zimmer

Edouard Bertrand

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

Florian Mueller
email: fmueller@pasteur.fr

FISH-quant v2: a scalable and modular tool for smFISH image analysis

Keywords: image analysis, RNA localization, transcription, smFISH

INTRODUCTION

Regulation of gene expression is essential for a cell to fulfill its basic functions, and its dysregulation can lead to serious failures at the cellular, tissular and organism level. Transcription levels are not only tightly regulated, but for many genes it has now been demonstrated that their transcripts accumulate in specific regions in the cell, thereby producing intricate localization patterns. Such subcellular targeting of mRNAs is thought to play an important role for the spatial control of gene expression and improper RNA trafficking is linked to an increasing number of diseases [START_REF] Buxbaum | In the right place at the right time: visualizing and understanding mRNA localization[END_REF][START_REF] Chin | RNA localization: making its way to the center stage[END_REF]. However, the function and mechanisms of RNA localization are not fully understood and we still lack a view of this process at the transcriptomic scale.

RNA abundance and localization can be studied at a large scale by image-based assays, where individual mRNA molecules are visualized by single-molecule Fluorescence in situ hybridization (smFISH). This technique allows for the detection of individual mRNA molecules in their native cellular environment [START_REF] Raj | Imaging individual mRNA molecules using multiple singly labeled probes[END_REF][START_REF] Tsanov | SmiFISH and FISH-quant: a flexible single RNA detection approach with super-resolution capability[END_REF]) by targeting each mRNA with several fluorescently labeled oligonucleotides. Many variants of this method exist, with optimizations regarding signal-to-noise ratio (SNR), experimental protocol, targeting specificity, scalability, automatization, and cost (for review, see [START_REF] Pichon | A growing toolbox to image gene expression in single cells: sensitive approaches for demanding challenges[END_REF]. Furthermore, an increasing number of multiplexing methods have also been proposed over the last years, enabling the simultaneous imaging of up to 10,000 RNA species in cells and tissues [START_REF] Moffitt | RNA imaging with multiplexed error-robust fluorescence in situ hybridization (MERFISH)[END_REF]Eng et al. 2019). Usually, smFISH experiments are complemented by the use of one or several fluorescent markers highlighting relevant compartments in the cell, such as the nucleus, the cytoplasm or any organelle that might serve as a reference, depending on the focus of the study.

These scalable imaging techniques produce extremely large and complex image data sets exploring spatial distributions of large portions of the transcriptome. While largescale imaging methods provide a systematic tool to understand RNA localization at a systems level, they come at a price: the need for fully automated, robust image analysis and user-friendly software tools to analyze such data sets and to fully exploit their potential [START_REF] Pichon | A growing toolbox to image gene expression in single cells: sensitive approaches for demanding challenges[END_REF][START_REF] Das | Intracellular mRNA transport and localized translation[END_REF].

Several specifications can be defined a priori for such an analysis tool. It should be simple enough to be mastered by non-experts, especially noncoders. Yet, it should be flexible enough to address different experimental designs and rely on a common algorithmic backbone. With the same modules, users should be able to both perform a high content screening analysis in a remote cluster, and a local analysis of a single image. Finally, the software should integrate the latest generation of computer vision algorithms, in particular deep-learning-based methods for image segmentation [START_REF] Ronneberger | U-Net: convolutional networks for biomedical image segmentation[END_REF][START_REF] Falk | U-Net: deep learning for cell counting, detection, and morphometry[END_REF][START_REF] Stringer | Cellpose: a generalist algorithm for cellular segmentation[END_REF].

Here, we introduce a Python-based version of our widely adopted software package FISH-quant [START_REF] Mueller | FISH-quant: automatic counting of transcripts in 3D FISH images[END_REF] for the analysis of smFISH images. Contrary to the first version of FISH-quant in Matlab, we address and improve on each of the specifications mentioned above. The switch to Python allows us to develop a flexible, free and fully open-source software. FISH-quant v2 enjoys a better integration to other open source tools and frameworks, from data analysis to web-based user interaction. Importantly, FISH-quant v2 facilitates the use of machine learning or deep learning algorithms with the import of dedicated packages, such as scikit-learn [START_REF] Pedregosa | Scikitlearn: machine learning in python[END_REF] or TensorFlow [START_REF] Abadi | TensorFlow: a system for large-scale machine learning[END_REF]. We also improve the scalability and the modularity of the package: the software has now been applied to several High Content Screening projects [START_REF] Chouaib | A dual protein-mRNA localization screen reveals compartmentalized translation and widespread co-translational RNA targeting[END_REF][START_REF] Pichon | The kinesin KIF1C transports APC-dependent mRNAs to cell protrusions[END_REF][START_REF] Safieddine | A choreography of centrosomal mRNAs reveals a conserved localization mechanism involving active polysome transport[END_REF]. Lastly, by using ImJoy [START_REF] Ouyang | ImJoy: an open-source computational platform for the deep learning era[END_REF]), a recently developed data analysis framework, we provide web-based graphical user interfaces (GUI) for both launching image analysis and downstream analysis of the results, and the computation can be performed locally or seamlessly scale to powerful remote computing servers.

RESULTS

The analysis of smFISH images aims at localizing and counting individual RNAs with respect to single cells and other subcellular landmarks. It typically encompasses a sequence of interconnected steps: (i) segmenting cells and the relevant cellular compartments such as nuclei (depending on the focus of the study and the markers used), (ii) detecting isolated and clustered RNA molecules, (iii) assignment of spots to cells, and (iv) analysis of expression levels and RNA localization patterns [START_REF] Battich | Image-based transcriptomics in thousands of single human cells at single-molecule resolution[END_REF][START_REF] Mueller | FISH-quant: automatic counting of transcripts in 3D FISH images[END_REF][START_REF] Stoeger | Computer vision for image-based transcriptomics[END_REF][START_REF] Tsanov | SmiFISH and FISH-quant: a flexible single RNA detection approach with super-resolution capability[END_REF][START_REF] Samacoits | A computational framework to study sub-cellular RNA localization[END_REF], potentially in combination with other phenotypic features [START_REF] Battich | Control of transcript variability in single mammalian cells[END_REF][START_REF] Safieddine | A choreography of centrosomal mRNAs reveals a conserved localization mechanism involving active polysome transport[END_REF].

Overview of existing analysis solutions

While several tools exist for each of these steps, there is currently-to our knowledge-no tool available that permits performing the entire analysis in one framework (see Supplemental Note 4). A complete analysis pipeline has then to be built by mixing these tools and requires some in-house developments, which can be daunting for nonspecialists and may provide solutions that are unstable and difficult to scale.

For the first step of object segmentation, deep-learning has become the method of choice with dramatic improvements in segmentation accuracy as compared to traditional methods. Several approaches exist that allow segmentation of cells and/or nuclei with minimal adjustment on new data sets, thanks to optimized models and large and diverse training data [START_REF] Schmidt | Cell detection with star-convex polygons[END_REF][START_REF] Hollandi | nucleAIzer: a parameterfree deep learning framework for nucleus segmentation using image style transfer[END_REF][START_REF] Lalit | Embedding-based instance segmentation in microscopy[END_REF][START_REF] Stringer | Cellpose: a generalist algorithm for cellular segmentation[END_REF]. The second step, fluorescence spot detection, has been addressed by a number of approaches in the literature, and more recently solutions specifically adapted to smFISH have been proposed. RS-FISH allows robust and accurate detection of fluorescent spots in 2D and 3D through radial symmetry but requires parameter tuning before being scaled to a large set of images [START_REF] Bahry | RS-FISH: precise, interactive and scalable smFISH spot detection using radial symmetry[END_REF]. DeepLink is a parameter-free deep-learning-based method, but is currently only available for 2D data and might require retraining [START_REF] Eichenberger | deepBlink: threshold-independent detection and localization of diffraction-limited spots[END_REF]. Lastly, assigning spot counts to segmentation results and the subsequent analysis of RNA levels and/or RNA localization requires custom-written code [START_REF] Stoeger | Computer vision for image-based transcriptomics[END_REF][START_REF] Samacoits | A computational framework to study sub-cellular RNA localization[END_REF].

General image analysis tools such as CellProfiler [START_REF] Mcquin | CellProfiler 3.0: next-generation image processing for biology[END_REF] permit us to establish an analysis framework daisy-chaining some of these analysis steps, but do not permit us to perform the entire analysis. A number of approaches, specifically dedicated to the analysis of smFISH are available. In our own software FISH-quant v1 [START_REF] Mueller | FISH-quant: automatic counting of transcripts in 3D FISH images[END_REF]) and also [START_REF] Stoeger | Computer vision for image-based transcriptomics[END_REF], the core of the analysis was performed in Matlab while cell segmentation was performed with the Python-based CellProfiler. DypFISH [START_REF] Savulescu | Interrogating RNA and protein spatial subcellular distribution in smFISH data with DypFISH[END_REF]) permits the study of the spatial distribution of mRNAs and proteins of micropatterned cells, mixing tools implemented in Python and Icy [START_REF] De Chaumont | Icy: an open bioimage informatics platform for extended reproducible research[END_REF]. Lastly, StarFISH [START_REF] Perkel | Starfish enterprise: finding RNA patterns in single cells[END_REF]) is an ongoing software FISH-quant v2 www.rnajournal.org 787 development mainly aiming at solving problems related to multiplex smFISH data for application in spatial transcriptomics.

FISH-quant v2: a complete toolbox for smFISH analysis

While an impressive range of methods already exists, a unified framework is lacking, which prevents users, especially non-specialist, from performing their smFISH analysis. To address this, we designed FISH-quant v2 to fulfill the above-described requirements in a flexible and efficient way. This version is entirely open-source and hosted on GitHub under the FISH-quant organization (Fig. 1, https ://github.com/FISH-quant). Using a GitHub organization allowed us to provide dedicated repositories with well defined and dedicated scope. Further, it gives the flexibility for future extension where new projects can be integrated as new, independent repositories, without affecting and complexifying the already existing code. The user can choose the adequate code for the analysis needs, without the overhead of installing unnecessary packages.

This GitHub organization is organized in several resources with dedicated repositories and documentation. First We recently showed how such images can be used to develop and validate analysis pipelines with the goal to quantify such intracellular RNA distributions [START_REF] Samacoits | A computational framework to study sub-cellular RNA localization[END_REF][START_REF] Dubois | A deep learning approach to identify mRNA localization patterns[END_REF]. Fourth, ImJoy plugins [START_REF] Ouyang | ImJoy: an open-source computational platform for the deep learning era[END_REF]) provide a graphical user-interface for the most commonly used workflows, and an interactive tutorial that can also run directly without local installation. Lastly, code from future projects either using or further improving FISH-quant will also be hosted here, creating a valuable, centralized resource for the community. A landing page (https://fish-quant.github.io/) directs new users to the most relevant resource for their analysis needs.

Big-FISH: Python package for smFISH analysis

We chose Python for the implementation of the core analysis package for several reasons: it allows the development of a free and fully open-source software, it provides established libraries for data and image analysis and is the language of choice for deep-learning implementations. Lastly, it can be interfaced with other tools and frameworks, from data analysis to web design, for instance with ImJoy [START_REF] Ouyang | ImJoy: an open-source computational platform for the deep learning era[END_REF] to provide interactive tools for user interaction and data inspection.

Our Python package includes several independent subpackages fitting the described workflow (see Materials and Methods for more details): preprocessing, segmentation, detection, and analysis. We designed each subpackage with clearly defined input and output data formats, which will be automatically checked. This then allows using each of these packages independently in a modular fashion. Users can thus create a customized analysis workflow, starting from preprocessing of images to statistical interpretation of results. These workflows can be implemented in Python and Bash scripts and run both on local and remote computational resources. The modular design also permits the easy integration of external methods, for instance, a new segmentation method can be combined with our spot detection algorithm. Lastly, we provide a subpackage to visualize the results of each intermediate step in the analysis workflow and thus provide valuable visual quality control.

Here, we will only provide an overview of these subpackages (Fig. 2). For a more detailed description of algorithms and methods, we refer to the documentation (https://bigfish.readthedocs.io/en/stable/) and the dedicated tutorials (https://github.com/fish-quant/big-fish-examples). These tutorials can be run directly in the browser with provided test data, and thus allow new users to immediately test these tools. The described methods were developed and validated with the data from two large-screen smFISH studies [START_REF] Chouaib | A dual protein-mRNA localization screen reveals compartmentalized translation and widespread co-translational RNA targeting[END_REF][START_REF] Safieddine | A choreography of centrosomal mRNAs reveals a conserved localization mechanism involving active polysome transport[END_REF]) (see Materials and Methods).

For image handling and preprocessing, we implemented a number of different utility functions to read, write, normalize, cast, filter, and project images. Different image file formats are natively supported and both 2D and 3D images can be processed.

The detection subpackage implements the methods required to detect spots in 2D or 3D images (Figs. 2, 3A-E). An important aspect of the detection subpackage is its ability to detect spots without setting any pixel intensity threshold. We implemented a method to automatically infer this threshold from the image. The curve describing the number of detected spots as a function of the intensity thresh-old (Fig. 3A,B) has an elbow shape, resulting from the superposition of the fast decreasing false positive detections (low intensity noise) and the slowly decreasing true positives. The threshold selected corresponds to the kink in the elbow, and corresponds thus to the highest threshold outside the high-noise regime. In order to validate this approach, we simulated realistic smFISH images with varying noise levels (Fig. 3A,B; Supplemental Note 1). We found that our method only leads to a moderate over-estimation of detected spots (<5%-10%) for images with moderate to high SNR values (>5). Such automatization overcomes human intervention and allows scaling to large data sets, such that the subpackage can process thousands of images. While initially designed to detect individual mRNAs, the same methods can also be used to detect other spot-like structures [START_REF] Safieddine | A choreography of centrosomal mRNAs reveals a conserved localization mechanism involving active polysome transport[END_REF]), such as centrosomes, Pbodies, etc (Fig. 3E).This subpackage further permits us to perform localization of RNAs with subpixel accuracy by using a Gaussian fitting [START_REF] Mueller | FISH-quant: automatic counting of transcripts in 3D FISH images[END_REF]. Lastly, we provide the possibility to perform a colocalization analysis between spot detection performed in multiple channels [START_REF] Cornes | piRNAs initiate transcriptional silencing of spermatogenic genes during C. elegans germline development[END_REF].

Strong local accumulation of RNAs, for example, active transcription sites, RNA foci, or areas of local translation [START_REF] Chouaib | A dual protein-mRNA localization screen reveals compartmentalized translation and widespread co-translational RNA targeting[END_REF], can lead to an underdetection since such accumulations are counted as single RNAs. For such cases, we provide tools to decompose these dense regions and estimate the number of spots based on our earlier work (Fig. 3C; [START_REF] Samacoits | A computational framework to study sub-cellular RNA localization[END_REF]. We validated this approach again on simulated data (Fig. 3D; see Supplemental Note 1), and found consistent performance across relevant noise levels.

The segmentation subpackage contains several algorithms and utility functions for segmentation and post-processing. It provides deep-learning-based approaches to segment cells and nuclei (Figs. 2, 3F,G; Supplemental Note 2). Furthermore, we provide postprocessing tools to refine and clean the segmentation result, such as boundary smoothing, removal of small objects or filling of small holes. Lastly, morphological properties, such as the area of cells, nuclei or protrusions, can be computed for these components (Supplemental Note 3).

The cell matching subpackage allows combining results from detection and segmentation, permitting us to analyze RNA abundance and distribution at the single-cell level. Detected spots can be assigned to a specific region of interest, for instance, a cell or a nucleus. Using the same method, RNA clusters can be assigned to a nucleus and thus be considered as transcription sites. RNA expression levels are extracted within this subpackage, as this is usual-ly the minimum information that is extracted from this kind of image.

The localization feature extraction subpackage permits the extraction of further information to study the subcellular spatial distribution of mRNA molecules. It gathers methods to format spot positions and coordinates of cellular landmarks and compute several spatial features at the singlecell level (Fig. 2; Supplemental Note 3). These features allow a statistical description of the cell population (Pichon et

Sim-FISH: simulation of smFISH images and RNA localization patterns

Simulations can be used to validate different steps of the analysis pipeline, ranging from the spot detection [START_REF] Tsanov | SmiFISH and FISH-quant: a flexible single RNA detection approach with super-resolution capability[END_REF]) to a statistical framework to quantitatively study RNA localization [START_REF] Samacoits | A computational framework to study sub-cellular RNA localization[END_REF][START_REF] Dubois | A deep learning approach to identify mRNA localization patterns[END_REF]). As mentioned above, we validated both our spot detection and our decomposition method for dense regions with this package (Fig. 3A-D; Supplemental Note 1). We simulate realistic smFISH images in three steps (see Supplemental Note 1). First, we randomly generate 2D or 3D spot coordinates, which can be random or display a specific subcellular RNA localization pattern. Further, clustered RNAs can be added. Second, we simulate a realistic image from these coordinates by modeling a RNA spot with a Gaussian function. Third, we add a noisy background to this image.

ImJoy: interactive user interfaces and data exploration

Our Python core analysis package provides flexibility and scalability since its components can be adapted to the specific analysis need of a given project. However, they require at least a minimum knowledge of Python to establish a complete workflow by using the provided tutorials.

To provide simpler access for users with no computational background and no programming skills, we implemented several plugins with graphical user interfaces for our computational platform ImJoy [START_REF] Ouyang | ImJoy: an open-source computational platform for the deep learning era[END_REF]). These plugins provide the most commonly used analysis workflow, as we determined from the usage of the Matlab version of FISH-quant, and will thus be suited for a large number of use cases (Fig. 4). First, a plugin to perform deep-learning-based segmentation. This is currently built on top of CellPose [START_REF] Stringer | Cellpose: a generalist algorithm for cellular segmentation[END_REF]), but thanks to our modular design, this can be easily exchanged if more performant methods are available in the future. Second, detection of both isolated and clustered RNA. Detection results can be conveniently inspected with the Kaibu image viewer plugin in ImJoy and different detection settings interactively investigated. Batch processing of entire folders is also possible. Lastly, detection results can be assigned to segmented cells and nuclei. We provide an interactive demo version of this plugin that can run directly in the browser without any local installation (https://fishquant.github.io/fq-interactive-docs/#/fq-imjoy).

Using ImJoy provides several advantages beyond simply providing a user interface. Due to its distributed design that separates GUI from computation plugins, it natively supports user-friendly remote computing which allows access to massive data storage and powerful computation resources including GPUs. ImJoy is a browser-based app where the user-interface plugin is implemented with JavaScript/CSS/HTML. ImJoy then transparently calls the computation functions in the Big-FISH package running on a Python plugin engine (e.g., Jupyter server) to perform the actual smFISH analysis task (Fig. 4). While this plugin can run on a local workstation, it can be executed on a computational cluster or even in the cloud or seamlessly switching between them. This is illustrated by the demo version, where the engine is running on Binder (Project [START_REF] Jupyter | Binder 2.0-reproducible, interactive, sharable environments for science at scale[END_REF]. Once the plugin engine is installed on the remote resource, the end-user can connect with ImJoy and will be confronted with the same interface, independently of where the analysis is actually performed. Interestingly, this front-end interface can also be opened with mobile devices, providing easy access.

ImJoy plugins implemented in JavaScript not only provide modern and reactive user-interfaces, but also profit from the extensive JavaScript data visualization libraries to build interactive data-inspection tools. Such interactivity is becoming increasingly important, especially when large and complex data sets are analyzed where static plots are too limited. As a case example, we provide an interactive t-SNE plot for the data shown in Figure 2 (https://fish-quant.github.io/fq-interactive-docs/ #/rnaloc-tsne). This plugin can be run without local installation and enables the user to explore and interact with these complex data.

Case studies

We developed and validated FISH-quant v2 for two largescale smFISH studies (see Materials and Methods). These two examples are typical use cases and exemplify the kind of quantitative results provided by this software. In [START_REF] Chouaib | A dual protein-mRNA localization screen reveals compartmentalized translation and widespread co-translational RNA targeting[END_REF], we performed a high-content screen in HeLa cells and analyzed 10,000 segmented cells. FISH-quant v2 was used for spot detection, cell segmentation and the computation of localization features that allowed us to apply supervised and unsupervised machine learning to identify localization patterns and classify single cells into predefined pattern classes. We observed several distinct mRNA localization patterns, including RNA accumulating (i) in foci, (ii) in cytoplasmic protrusions, (iii) in the perinuclear area (which could be subdivided in endosomal, RE, Golgi and centrosome associate), (iv) forming a rim at the nuclear edge, or (v) inside the nucleus (Fig. 2). Interestingly, automated classification done on a single-cell level revealed a high degree of cell-to-cell heterogeneity in RNA localization, with 10% to 80% of the cells displaying the expected pattern depending on the RNA (Fig. 5A). In addition, for each pattern, only a fraction of the mRNA appeared to localize, revealing a high degree of plasticity in RNA localization mechanisms. This appears to be specific to cell lines as RNA localization in embryos is usually much more stereotyped. We also quantified how translation inhibition affected RNA localization and found that most mRNAs localize in a translation-dependent manner, which is unexpected (Fig. 5B). This also enabled us to discover translation factories, small cytoplasmic structures where specific mRNAs accumulate to be translated.

In [START_REF] Safieddine | A choreography of centrosomal mRNAs reveals a conserved localization mechanism involving active polysome transport[END_REF], we studied RNA localization at centrosomes (3600 images and 54,000 cells). Here, we added an automated detection for centrosomes (Fig. 3E), and implemented localization features describing this localization pattern (Fig. 5C). This enabled us to discover a family of eight centrosomal mRNAs whose localization to centrosome is cell cycle dependent and conserved from humans to drosophila.

Altogether, these analyses demonstrate the power of FISH-quant v2 in processing large smFISH data sets, and classifying RNA localization patterns in an automated way.

DISCUSSION

Here, we present FISH-quant v2, a user-friendly Pythonbased software for the complete analysis of smFISH images. It is built around a core-analysis package, implemented following rigorous software development guidelines, with detailed interactive documentation and tutorials. This package consists of several interchangeable modules permitting the construction of highly flexible workflows for specific analysis needs. For standard workflows, we provide user interfaces in ImJoy accessible to biologists without programming skills, which can be used locally or scaled to larger remote computational resources. Finally, FISHquant hosts a simulation package to generate smFISH images with nonrandom intracellular RNA localization patterns. These simulated images can be used to develop and evaluate analysis pipelines to study such RNA localization [START_REF] Samacoits | A computational framework to study sub-cellular RNA localization[END_REF][START_REF] Dubois | A deep learning approach to identify mRNA localization patterns[END_REF]. As demonstrated in two recently published studies [START_REF] Chouaib | A dual protein-mRNA localization screen reveals compartmentalized translation and widespread co-translational RNA targeting[END_REF][START_REF] Safieddine | A choreography of centrosomal mRNAs reveals a conserved localization mechanism involving active polysome transport[END_REF], FISH-quant v2 can be used for large screening data sets thanks to its scalability. Spot detection, segmentation, feature extraction and pattern recognition can be performed over thousands of cells without fine-tuning parameters for every image.

We designed FISH-quant v2 based on the successful previous implementation in Matlab [START_REF] Mueller | FISH-quant: automatic counting of transcripts in 3D FISH images[END_REF] integrating new features and user feedback we obtained from several projects over several years. The entire core package is written in Python since this allowed us to address the above-mentioned requirements for a smFISH analysis tool. We use established scientific libraries (see Materials and Methods), and keep these dependencies to a minimum facilitating installation, maintenance and the integration with other analysis frameworks. These libraries are developed, validated and maintained by a large scientific community, ensuring long-term support and availability. We further use strict version control, guaranteeing reproducibility. Lastly, all dependencies, as well as FISH-quant v2, are open-source, thus can be used free of charge, both on local and remote computational infrastructures, and thus analysis can easily be scaled to larger data volumes.

The organization of the analysis subpackages in the core package matches key steps in smFISH image analysis, with a special focus on flexibility. All steps (preprocessing, RNA detection, segmentation as well as data inspection and analysis) can be run independently or replaced by external code, by respecting a strict data format. This allows FISHquant to be adapted to the respective analysis needs, and build custom workflows.

While this flexibility is important, many users require a standard workflow and do not have programming experience. For these cases, we provide ImJoy plugins with a convenient user interface running in the browser [START_REF] Ouyang | ImJoy: an open-source computational platform for the deep learning era[END_REF]). These interfaces are built with modern web libraries and are thus intuitive, and no experience in Python is required to analyze data. Lastly, these ImJoy plugins can be readily extended by more experienced users to further adapt them to their needs. A detailed documentation and interactive tutorial further help new users to get started quickly.

In summary, we present with FISH-quant v2 a rigorously validated analysis platform for smFISH data, developed to match the analysis requirements of large data sets. Its modularity permits the creation of flexible workflows ranging from the analysis of small data sets with the help of a graphical user-interface to custom-tailored investigation of large-scale screens requiring computational clusters.

MATERIALS AND METHODS

Python core packages

The repository Big-FISH contains the Python code used for the actual analysis. It is organized in several subpackages performing dedicated steps:

-I/O operations, images preprocessing and (bigfish.stack) -mRNA spot detection (bigfish.detection) -nucleus and cell segmentation (bigfish.segmentation) -post-processing and analysis of results from different channels, such as the merging of RNA detections and segmentation masks or colocalization analysis (bigfish.multistack)

-feature computation, point cloud analysis and classification (bigfish.classification)

-visual reports of the obtained results (bigfish.plot)

-application of deep learning algorithms for segmentation (bigfish.deep_learning)

The repository Sim-FISH contains the Python code used for simulations. It includes several modules to generate 3D spots coordinates (both random and with a specific subcellular localization pattern). From these coordinates, simulated smFISH images with a noisy background can be generated. Dependencies are limited to standard Python scientific libraries: scientific computing (numpy [START_REF] Harris | Array programming with NumPy[END_REF] and SciPy [START_REF] Virtanen | SciPy 1.0: fundamental algorithms for scientific computing in Python[END_REF]), data wrangling (pandas [McKinney 2010]), image analysis (scikit-image [van der Walt et al. 2014]), visualization (matplotlib [START_REF] Hunter | Matplotlib: a 2D graphics environment[END_REF]]), parallel computing (joblib, https://github.com/joblib/joblib) and machine learning (scikit-learn [START_REF] Pedregosa | Scikitlearn: machine learning in python[END_REF], TensorFlow [START_REF] Abadi | TensorFlow: a system for large-scale machine learning[END_REF]).

The GitHub repositories are using continuous integration providing increased robustness of the released code, through unitary testing, version control and automatically generated up-to-date documentation. Packages are hosted under a BSD 3-Clause License.

Example data sets

Two data sets were used for the development and validation of FISH-quant. First, from a screen studying local translation and consisting of 526 fields of view (DAPI and smFISH channels) from 57 separate experiments (27 different mRNAs under different experimental conditions [START_REF] Chouaib | A dual protein-mRNA localization screen reveals compartmentalized translation and widespread co-translational RNA targeting[END_REF]). For this screen, 3D images with a z-spacing of 0.3 µm were acquired on two different systems: (i) a Zeiss AxioimagerZ1 wide-field microscope equipped with a motorized stage, a camera scMOS ZYLA 4.2 MP, using 63× and 100× oil objectives, (ii) Nikon Ti fluorescence microscope equipped with ORCA-Flash 4.0 digital camera (HAMAMATSU). Second, from a screen focusing on local translation of centrosomal mRNAs. The data set consisted of 3678 fields of view (Dapi, smFISH, CellMask and GFP channels) from 218 experiments [START_REF] Safieddine | A choreography of centrosomal mRNAs reveals a conserved localization mechanism involving active polysome transport[END_REF]). 3D images were acquired with an automated spinning disk microscope (Opera, PerkinElmer), equipped with a 63× water objective. Z-spacing was 0.3 µm.

FISH-quant v2

www.rnajournal.org 793

  , a Python package (Big-FISH) providing the core code for performing computation and analysis. Second, detailed interactive examples with test data for each analysis step implemented in Jupyter notebooks. These examples can be run directly on Binder (Project Jupyter et al. 2018), a free and reproducible Jupyter notebook service, without local installation. Third, a repository containing code to simulate different subcellular RNA localization patterns (Sim-FISH).

FIGURE 1 .

 1 FIGURE 1. Organization of FISH-quant. FISH-quant is hosted on GitHub and consists of several interconnected repositories. The Python core package contains the entire analysis code, which is used by both the ImJoy plugins and the example and tutorial repository.

FIGURE 2 .

 2 FIGURE 2. Big-FISH: the core analysis Python analysis package. (Upper part) Main modules illustrated with a typical analysis workflow. Shown are also the inputs and outputs that are created at the different steps. (Lower part) As a final result of the analysis of Big-FISH, each cell is described with a set of features reflecting RNA abundance and localization. These features can then be used to perform analysis on the cell population. Shown are results from our RNA localization screen where cells are grouped based on their RNA localization pattern (Chouaib et al. 2020). The t-SNE plot projects 15 localization features for smFISH experiments against 27 different genes. Each dot is one cell. The color-coded dots are manual annotations of six different localization patterns. Images are examples of individual cells displaying a typical localization pattern of this region of the t-SNE plot.

  FIGURE 3. (A) Automated spot detection. Simulated image (left) and detection results (right) with detected spots in red and ground truth in white. (B) Elbow curve used for automated threshold setting, red dot indicates identified intensity threshold. (C) Decomposition of dense regions. Simulated image (left) and decomposition results (right) with detected spots in red and ground truth in white. Number of simulated and detected spots are shown in white and red, respectively. (D) Algorithm to decompose dense regions was evaluated with 100 simulated images containing a cluster of 15 spots and different noise levels. (E) Example of automated detection of BICD2 mRNAs (left) and centrosome (right) in HeLa cells. (F) Example of nucleus segmentation from a DAPI image. (G) Example of cell segmentation from a CellMask image.

FIGURE 4 .

 4 FIGURE 4. ImJoy. Schematic view of Imjoy's architecture. ImJoy's core is a Progressive Web App whose functionalities are provided by plugins that can be written in different programming languages. ImJoy can perform computations in the browser (including offline), locally or remotely via plugin engines.

  FIGURE 5. (A) Heatmap depicting the fraction of cells classified in the indicated pattern, for the different genes analyzed by the automated pipeline. (B) Impact of treatment with translational inhibitor puromycin on the number of detected RNA clusters. HMMR shows a similar number of clusters, while all other genes have significantly fewer, indicating an implication of translation in cluster formation. (C) Proportion of mRNAs within 2000 nm of a centrosome. Distance threshold was empirically defined as the typical distance between clustered RNAs and the centrosomes. Compared are untreated cells, and cells treated with two different translation inhibitors: cycloheximide, blocking ribosome elongation, or puromycin, inducing premature chain termination. BICD2 has a centrosomal localization pattern, while TRIM59 is a negative control with a random intracellular localization. Results are displayed with different treatments.

  Cold Spring Harbor Laboratory Press on January 17, 2023 -Published by rnajournal.cshlp.org Downloaded from

Cold Spring Harbor Laboratory Press on January 17, 2023 -Published by rnajournal.cshlp.org Downloaded from

RNA (2022) Vol. 28, No. 6 Cold Spring Harbor Laboratory Press on January 17, 2023 -Published by rnajournal.cshlp.org Downloaded from

www.rnajournal.org 789 Cold Spring Harbor Laboratory Press on January 17, 2023 -Published by rnajournal.cshlp.org Downloaded from

www.rnajournal.org 791 Cold Spring Harbor Laboratory Press on January 17, 2023 -Published by rnajournal.cshlp.org Downloaded from

ACKNOWLEDGMENTS

This work was funded by the ANR (ANR-19-CE12-0007) and Institut Pasteur and 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). Furthermore, we also acknowledge France-BioImaging infrastructure supported by the French National Research Agency (ANR-10-INBS-04).

Received December 2, 2021; accepted February 19, 2022.

DATA DEPOSITION

The entire code for the analysis described in this paper is available on GitHub: https://github.com/fish-quant. This study includes no data deposited in external repositories.

SUPPLEMENTAL MATERIAL

Supplemental material is available for this article.


Peter Naylor

Tristan Lazard

Guillaume Bataillon

Marick Laé

Anne Vincent-Salomon

Anne-Sophie Hamy

Fabien Reyal

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

Elaine Wan

Ling Chan

Jónathan Heras

Prediction of Treatment Response in Triple Negative Breast Cancer From Whole Slide Images

Keywords: breast cancer, digital pathology, whole slide images, treatment response, cross-validation, deep learning, multiple instance learning, small sample size

The automatic analysis of stained histological sections is becoming increasingly popular. Deep Learning is today the method of choice for the computational analysis of such data, and has shown spectacular results for large datasets for a large variety of cancer types and prediction tasks. On the other hand, many scientific questions relate to small, highly specific cohorts. Such cohorts pose serious challenges for Deep Learning, typically trained on large datasets. In this article, we propose a modification of the standard nested crossvalidation procedure for hyperparameter tuning and model selection, dedicated to the analysis of small cohorts. We also propose a new architecture for the particularly challenging question of treatment prediction, and apply this workflow to the prediction of response to neoadjuvant chemotherapy for Triple Negative Breast Cancer.

INTRODUCTION

Context

Breast cancer (BC) is the most common cancer in women and the leading cause of cancer deaths with 18.2% of deaths among female cancer patients and 8% among all cancer patients (Institut National Du Cancer, 2019). Out of the four main breast cancer types, Triple Negative Breast Cancer (TNBC) represents 10% of all BC patients. This group has the worst prognostic with a five-year survival rate of around 77 percent versus 93 percent for the others. Currently no specialised treatments exists and the standard procedure consists in administrating neoadjuvant or adjuvant chemotherapy [START_REF] Sakuma | Pathological Tumor Response to Neoadjuvant Chemotherapy Using Anthracycline and Taxanes in Patients with Triple-Negative Breast Cancer[END_REF]. TNBC research is still a very active field of study [START_REF] Foulkes | Triple-negative Breast Cancer[END_REF] and on the one hand, most works have focused on stratifying cohorts based on molecular and biological profiles [START_REF] Lehmann | Refinement of Triple-Negative Breast Cancer Molecular Subtypes: Implications for Neoadjuvant Chemotherapy Selection[END_REF]. We, on the other hand, tackle the problem of predicting the response variable in a TNBC neoadjuvant chemotherapy (NACT) cohort from a histological needle-core biopsy section from the primary tumor prior to treatment. In contrast to most of the effort in cancer research, which is driven by the analysis of sequencing data, our study is based solely on the histological image data prior to treatment.

Each histological sample corresponds to tissue slides encompassing the tumor and its surrounding, stained with agents in order to highlight specific structures, such as cells, cell nuclei or collagen. The morphological properties of these elements and their spatial organisation have been linked to cancer subtypes, grades and prognosis. Even if pathologists have been trained to understand and report the evidence found in this type of data, the complexity, size and heterogeneity found in histological specimens make it highly unlikely that all relevant and exploitable patterns are known of today. Tissue images are informative about morphological and spatial patterns and are therefore inherently complementary to omics data.

Two major technological advances have triggered the emergence of the field of Computational Pathology: first, the arrival of new and powerful scanners replaced to some extent the use of conventional microscopes. Today, slides are scanned, stored and can be accessed rapidly at any moment [START_REF] Huisman | Creation of a Fully Digital Pathology Slide Archive by High-Volume Tissue Slide Scanning[END_REF]. This in turn has lead to the generation of large datasets that can also be analyzed computationally. The second element was the rise of new computer vision techniques.

Indeed, while the analysis of tissue slides has been of interest to the Computer Vision community for many years [START_REF] Bartels | Machine Learning in Quantitative Histopathology[END_REF], it is the advent of deep learning that has truly impacted the field. The advent of deep learning has stemmed a wide number of projects and investments. For visual systems, it is the combination of very large annotated dataset [START_REF] Jia | ImageNet: A Large-Scale Hierarchical Image Database[END_REF], hardware improvement and convolutional neural networks (CNN) that led to humanlike capabilities. These outbreaks in performance have led to the creation of many annotated datasets and to the application of CNN's to many tasks. However, for biomedical imaging in particular the application of CNN is not always straightforward:

1) The price for generating and annotating large biomedical dataset limits the progress of big data in this domain [START_REF] Gurcan | Histopathological Image Analysis: A Review[END_REF]. 2) In histopathology in particular, each individual sample can be very large, one sample can be up to 60 GB uncompressed. This leads to multiple issues, again linked with the time and price needed for annotation, but also for the subsequent analysis where ad hoc methods have to be used as an entire image does not comfortably fit in RAM.

3) The nature of the data is inherently complex, each biological sample has its own individual patterns to be differentiated with relevant pathological evidence. For histopathology data, samples have a very large inter-slide, but also intra-slide variability that make the apparent signal harder to detect.

In addition, the level of detail can be an additional difficulty: the relevant image features may be very fine grained, such as mitotic events, or very large such as the size of relevant image regions (necrotic, tumorous) [START_REF] Janowczyk | Deep Learning for Digital Pathology Image Analysis: A Comprehensive Tutorial with Selected Use Cases[END_REF].

In this paper, we set out to predict the response to NACT in TNBC from Whole Slide Images (WSI). Prediction of treatment response is one of the most difficult tasks in digital pathology [START_REF] Echle | Deep Learning in Cancer Pathology: a New Generation of Clinical Biomarkers[END_REF]. Unlike tasks like tumor detection and subtyping, for which high accuracies have been achieved in the past [START_REF] Echle | Deep Learning in Cancer Pathology: a New Generation of Clinical Biomarkers[END_REF], there are no known morphological phenotypes that would allow for the prediction of treatment response. Moreover, datasets with treatment response are usually much smaller than datasets acquired in clinical routine for tasks such as automatic grading, metastasis detection or subtype prediction.

The contributions of this article are:

1) We benchmark several state-of-the-art architectures with respect to their performance in treatment response prediction. 2) We propose a new architecture for the prediction of treatment response and demonstrate its efficiency. 3) We propose a suitable model selection procedure that can cope with small datasets and avoids re-training, which is particularly important for deep learning. We prove its validity and efficiency on simulated data.

The paper is organised as follows: in the next Section 1.2 we describe related work. In Section 2, we describe the methodological developments. In Section 2.1, we present the limits of the current validation procedures and our alternative method used in this study. We then introduce our histopathology dataset in Section 2.3. Section 2.4 is devoted to introducing the DNN architectures which will be applied to the TNBC cohort. In Section 3, we show our results on the simulated data for our validation procedure and the application of our DNN to the TNBC cohort. Finally in Section 4 we discuss our methods and results.

Related Work

Challenges in Computational Pathology

The field of research in computational pathology can be divided in three categories: 1) Preprocessing, in particular color normalisation which aims at reducing the bias introduced by staining protocols used in different centers [START_REF] Ruifrok | Quantification of Histochemical Staining by Color Deconvolution[END_REF][START_REF] Ehteshami Bejnordi | Stain Specific Standardization of Whole-Slide Histopathological Images[END_REF]. 2. Detection, segmentation and classification of objects of interest, such as regions [START_REF] Bejnordi | Deep Learning-Based Assessment of Tumor-Associated Stroma for Diagnosing Breast Cancer in Histopathology Images[END_REF][START_REF] Chan | Histosegnet: Semantic Segmentation of Histological Tissue Type in Whole Slide Images[END_REF][START_REF] Barmpoutis | Tertiary Lymphoid Structures (TLS) Identification and Density Assessment on H&E-stained Digital Slides of Lung Cancer[END_REF] and nuclei (Naylor et al., 2018;[START_REF] Graham | Hover-Net: Simultaneous Segmentation and Classification of Nuclei in Multi-Tissue Histology Images[END_REF].

3) The prediction of slide variables, such as presence and detection of disease [START_REF] Litjens | 1399 H&E-stained sentinel Lymph Node Sections of Breast Cancer Patients: the CAMELYON Dataset[END_REF][START_REF] Campanella | Clinical-grade Computational Pathology Using Weakly Supervised Deep Learning on Whole Slide Images[END_REF], grading [START_REF] Niazi | Visually Meaningful Histopathological Features for Automatic Grading of Prostate Cancer[END_REF][START_REF] Ström | Artificial Intelligence for Diagnosis and Grading of Prostate Cancer in Biopsies: a Population-Based, Diagnostic Study[END_REF], survival [START_REF] Zhu | WSISA: Making Survival Prediction from Whole Slide Histopathological Images[END_REF][START_REF] Courtiol | Deep Learning-Based Classification of Mesothelioma Improves Prediction of Patient Outcome[END_REF], gene expression [START_REF] Binder | Towards Computational Fluorescence Microscopy: Machine Learning-Based Integrated Prediction of Morphological and Molecular Tumor Profiles[END_REF][START_REF] Schmauch | A Deep Learning Model to Predict RNA-Seq Expression of Tumours from Whole Slide Images[END_REF], genetic mutations [START_REF] Coudray | Classification and Mutation Prediction from Non-small Cell Lung Cancer Histopathology Images Using Deep Learning[END_REF] or genetic signatures [START_REF] Kather | Pan-cancer Image-Based Detection of Clinically Actionable Genetic Alterations[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].

Of note, these different of Computer Vision tasks have varying degrees of difficulty. For instance, Deep Learning based tumor detection and subtyping can be achieved with high accuracies [AUC: 0.97-0.99 for detection, AUC: 0.85-0.98 for subtyping, [START_REF] Echle | Deep Learning in Cancer Pathology: a New Generation of Clinical Biomarkers[END_REF]]. This is not surprising, as these tasks rely on well-known visual cues. In contrast, prediction of treatment response is deemed one of the most difficult tasks in digital pathology [START_REF] Echle | Deep Learning in Cancer Pathology: a New Generation of Clinical Biomarkers[END_REF], because the related morphological patterns are unknown and could potentially be very complex. It is even not clear to which extent treatment response is actually predictable from WSI.

Pipelines for slide variable predictions are usually divided into several steps. Tiles are partitioned into smaller images, usually referred to as patches or tiles, which are then encoded by a DNN, often trained on ImageNet [START_REF] Zhu | WSISA: Making Survival Prediction from Whole Slide Histopathological Images[END_REF][START_REF] Courtiol | Deep Learning-Based Classification of Mesothelioma Improves Prediction of Patient Outcome[END_REF], as depicted in Figure 2. The training on ImageNet might be surprising at first sight, as the nature of the images are very different. In addition, ImageNet samples usually have a natural orientation, where the main object of interest is usually centered and scaled to fit in the image [START_REF] Raghu | NeurIPS 2019[END_REF]. In contrast, histopathology images have a rotationally invariant content with no prior regarding scale or positioning of the relevant structures. However, rotational invariance can be imposed (Naylor et al., 2019;[START_REF] Lafarge | Roto-Translation Equivariant Convolutional Networks: Application to Histopathology Image Analysis[END_REF], and in practice ImageNet based encodings are widely used and tend to perform very well.

After encoding of all tiles, each WSI is converted into a P × n i matrix where P is the encoding size and n i the number of tiles. The last step consists in aggregating tile level encodings to perform unsupervised or supervised predictions at the slide level [START_REF] Zhu | WSISA: Making Survival Prediction from Whole Slide Histopathological Images[END_REF][START_REF] Courtiol | Classification and Disease Localization in Histopathology Using Only Global Labels: A Weakly-Supervised Approach[END_REF][START_REF] Campanella | Clinical-grade Computational Pathology Using Weakly Supervised Deep Learning on Whole Slide Images[END_REF][START_REF] Courtiol | Deep Learning-Based Classification of Mesothelioma Improves Prediction of Patient Outcome[END_REF]Naylor et al., 2019).

Computational Pathology as a field has benefited from the generation of large annotated data sets, mostly with pixel-level annotations [START_REF] Kumar | A Dataset and a Technique for Generalized Nuclear Segmentation for Computational Pathology[END_REF][START_REF] Litjens | 1399 H&E-stained sentinel Lymph Node Sections of Breast Cancer Patients: the CAMELYON Dataset[END_REF]Naylor et al., 2018), or cell-level annotations [START_REF] Veta | Assessment of Algorithms for Mitosis Detection in Breast Cancer Histopathology Images[END_REF] for cell classification. The major resource for WSI with slide level annotations are the Cancer Genome Atlas (TCGA) and the Camelyon Challenge [START_REF] Litjens | 1399 H&E-stained sentinel Lymph Node Sections of Breast Cancer Patients: the CAMELYON Dataset[END_REF]. These public repositories are paralleled by many in-house datasets (thus not accessible to the public), some of which can be very large, namely in a screening context, e.g., [START_REF] Campanella | Clinical-grade Computational Pathology Using Weakly Supervised Deep Learning on Whole Slide Images[END_REF]. In most cases however, the datasets tend to be very small and fall therefore in the small n large p category. This is due to the fact that often the most interesting studies focus on particular molecularly defined cancer subtypes for which only small cohorts exist. In addition, collecting the output variable might be very challenging and time-consuming, if the project is not formulated in the context of Computer Aided Diagnosis. This is particularly true for treatment response prediction.

Challenges in Applying DNN to Small Datasets

For all supervised learning method it is custom to use a two step procedure for estimating the performance. After dividing your dataset into three categories: train, validation and test. The first step consists in performing model selection with the training and validation set. The second step simply involves evaluating the chosen model on the test set in order to assess an unbiased estimator of the performance [START_REF] Pereira | Machine Learning Classifiers and Fmri: a Tutorial Overview[END_REF]. This is however only possible if the three categories are large enough. When the validation set is too small and the discrepancy in the data too high, one could very easily over-fit or under-fit on the dataset [START_REF] Varoquaux | Cross-validation Failure: Small Sample Sizes lead to Large Error Bars[END_REF]. When the number of samples n is small, which is usually the case for biomedical data, alternative validation methods have to be found such as cross-validation (CV) and nested cross-validation (NCV). CV is mostly used for model selection or assessing performance. NCV is used when the model needs a tuning based on an external dataset, such as hyperparameter tunning for Support Vector Machines. Even if these methods have been debated [START_REF] Krstajic | Cross-validation Pitfalls when Selecting and Assessing Regression and Classification Models[END_REF][START_REF] Varoquaux | Cross-validation Failure: Small Sample Sizes lead to Large Error Bars[END_REF][START_REF] Wainer | Nested Cross-Validation when Selecting Classifiers Is Overzealous for Most Practical Applications[END_REF], they are widely accepted. These methods are explained in more details in Section 2.1.

Prediction of the Response to Neoadjuvant Chemotherapy in Triple Negative Breast Cancer

Neoadjuvant chemotherapy (NACT) responses varies among patients in TNBC and no clear biological signal has been shown. Survival in these cohorts have been correlated to the Residual Cancer Burden (RCB) variable [START_REF] Symmans | Measurement of Residual Breast Cancer burden to Predict Survival after Neoadjuvant Chemotherapy[END_REF] which can be used as a proxy for response. RCB is a pathological variable based on measurements of how much the primary tumor has shrunk and of the size of metastasis in axillary lymph nodes. Finding biological evidence to NACT response would allow for adequate and specific treatment, some histological variables have been found to be correlated with survival, such as the number Ki-67 positive cells [START_REF] Elnemr | Response of Triple Negative Breast Cancer to Neoadjuvant Chemotherapy: Correlation between Ki-67 Expression and Pathological Response[END_REF], tumor infiltrating lymphocytes [START_REF] Mao | The Prognostic Value of Tumor-Infiltrating Lymphocytes in Breast Cancer: a Systematic Review and Meta-Analysis[END_REF] and the Elston and Ellis grade [START_REF] Elston | Pathological Prognostic Factors in Breast Cancer. I. The Value of Histological Grade in Breast Cancer: Experience from a Large Study with Long-Term Follow-Up[END_REF]. Depending on the context, some alternative treatments have been found to help overall survival, such as those based on anthracycline and taxanes [START_REF] Sakuma | Pathological Tumor Response to Neoadjuvant Chemotherapy Using Anthracycline and Taxanes in Patients with Triple-Negative Breast Cancer[END_REF], carboplatin [START_REF] Pandy | Triple Negative Breast Cancer and Platinum-Based Systemic Treatment: a Meta-Analysis and Systematic Review[END_REF] or with olaparib and talazoparib [START_REF] Won | Triple-negative Breast Cancer Therapy: Current and Future Perspectives[END_REF]. Some treatments have emerged with targeted immunotherapy in combination with atezolizumab (anti-PD-L1 antibody) and nanoparticle albumin-bound (nab)paclitaxel [START_REF] Won | Triple-negative Breast Cancer Therapy: Current and Future Perspectives[END_REF]. Most of the studies for NACT responses have been performed in clinical practices and based on pathological variables [START_REF] Elnemr | Response of Triple Negative Breast Cancer to Neoadjuvant Chemotherapy: Correlation between Ki-67 Expression and Pathological Response[END_REF][START_REF] Gass | Prediction of Pathological Complete Response and Prognosis in Patients with Neoadjuvant Treatment for Triple-Negative Breast Cancer[END_REF][START_REF] Zhu | Predictors of Response and Survival Outcomes of Triple Negative Breast Cancer Receiving Neoadjuvant Chemotherapy[END_REF]. In addition, some studies have analysed sequencing and molecular profiles in order to better understand and stratify cohorts [START_REF] Lehmann | Refinement of Triple-Negative Breast Cancer Molecular Subtypes: Implications for Neoadjuvant Chemotherapy Selection[END_REF][START_REF] García-Vazquez | A Microrna Signature Associated with Pathological Complete Response to Novel Neoadjuvant Therapy Regimen in Triple-Negative Breast Cancer[END_REF][START_REF] Wang | Molecular Stratification within Triple-Negative Breast Cancer Subtypes[END_REF].

To the best of our knowledge, it remains unclear whether and to which extent NACT response can be predicted from biopsies taken prior to treatment, and only few works have addressed this question so far (Naylor et al., 2019;[START_REF] Ogier Du Terrail | Collaborative Federated Learning behind Hospitals' Firewalls for Predicting Histological Response to Neoadjuvant Chemotherapy in Triple-Negative Breast Cancer[END_REF].

MATERIALS AND METHODS

Validation Procedure

Here, we present our procedures to replace NCV in order to train DNN in a context of small n large p. We first explain crossvalidation (CV), nested cross-validation (NCV) and their limitations. We propose two different procedures, better suited and show their effectiveness on a small case study.

Cross-Validation

CV is a common procedure for model evaluation or model selection, specifically in situations where the data set is relatively small. CV divides the initial data set into k f folds, denoted F cv 1 , F cv 2 , . . . , F cv k f and runs algorithms on the data sets with one fold left out. We define, for all j, the set F cv -j , which is the union of all folds expect for fold j:

F cv -j ⋃ k∈[[1,k f ] ] k≠j F cv k 2.1.1.

Cross-Validation for Model Selection

CV can be used for model selection or model tuning. The procedure that returns a tuned model M will be notated f cv .

M f cv D ( ) (1)
We give the pseudo code in Algorithm 1, where H h 1 , . . . h i , ‥ is the set of hyperparameters (HP).

Algorithm 1. Model selection, f cv .

Cross-Validation for Model Evaluation

We can use CV to evaluate a given model and a HP set h. The procedure is similar to the pseudo code given in Algorithm 1, however we give in input a model and only one hyperparameter set and return t h . In this case, t h is an unbiased estimator of the performance, as no optimization of the hyperparameters took place. If however several sets of hyperparameters are tested as to minimize the accuracy measured by cross-validation, this accuracy is an overoptimistic estimation of the true accuracy. In order to get a realistic estimation of the accuracy, we therefore have to turn to nested cross-validation.

Nested Cross-Validation

NCV is a procedure that allows one to tune a model and effectively report an unbiased estimation of the performance of the tuned model.

Given sets of HP and a data set D, NCV corresponds to two nested loops of CV: The outer CV loop is for model evaluation, usually applied on test folds, sometimes referred to as outer folds. The inner CV loop is for model tuning, i.e., for each test fold, we perform a complete CV on the remaining data to correctly tune the model, and test the performance of the tuned model on data that has neither been used for training nor for HP tuning. We show the pseudo-code for NCV in Algorithm 2.

Algorithm 2. Nested cross-validation.

Another possible view is to see NCV as a simple CV for a model selection algorithm. For NCV, the model selection algorithm would be f cv .

It is important to note that as we are training DNN, we do not use a fixed hyperparameter set, H, but randomly generate the set as it has been shown that randomised search performs better [START_REF] Bergstra | Random Search for Hyper-Parameter Optimization[END_REF].

Limitations

DNN training suffers from inherent randomness, as the loss function is highly non-convex and possess many symmetries [START_REF] Bishop | Pattern Recognition[END_REF]. In addition, there are some stochastic differences between training runs, such as the random initialization of the weight parameters and the data shuffling, naturally leading to different solutions. Especially for small datasets, these stochastic variations lead to notable differences in performance when we repeat training with the same hyperparameters.

In the classical setting, CV provides us with a set of hyperparameters that lead to a model with optimal performance, as estimated in the inner loop. For DNN trained on small datasets, there is no guarantee that the same set of hyperparameters will lead to similar performance, and for this reason retraining is not guaranteed to lead to a very good solution.

Another problem with the retraining in line 10 in Algorithm 1 is the use of early stopping. Early stopping is a very powerful regularization procedures that choses experimentally the point between the under-and overfitting regime, but for this it requires a validation set. Early stopping would therefore not be applicable in the traditional CV-scheme with retraining.

Nested Ensemble Cross-Validation

Due to the incompatibility between NCV and early stopping we propose to modify the model selection procedure, i.e., function f cv shown in Algorithm 1. In particular we do not perform retraining and return an ensemble of the models used during CV. Similarly to NCV, we perform a CV where we propose to modify f cv into a better suited procedure, named f ecv (ensemble cross-validation), shown in Algorithm 3. The main difference between f cv and f ecv is that we remove the final model retraining, i.e., line 10 of Algorithm 1 and give back the full set of k f models trained for all folds for the maximizing hyperparameters; the prediction is obtained by ensembling these models.

The advantage of this procedure is that we omit the retraining step which allows us to use early stopping for all individual models. In addition, we add another level of regularization by the ensembling. Of note, f ecv can be used in an inner loop, too. This then leads to Nested Ensemble crossvalidation (NECV).

Nested Best Fold Ensemble Cross-Validation

Similarly to the procedure NECV, we can define another procedure that we name NBFCV. NBFCV is based on f bfcv

Simulations

In order to compare and validate our procedure presented in Section 2.1 we propose to conduct a series of simulation studies where the results will be given in Section 3.1. In particular, we wish to demonstrate that DNN training given a set of HP can lead to inconsistent models and that NCV therefore might provide under-performing models compared to NECV and NBFCV, the proposed validation procedures.

Data Set Simulation

We simulated a simple balanced binary data, of size N = 350 in a medium to high dimensional setting with p = 256. We have ∀i ∈ [(1, N)], Y i ∈ (-1, 1) and X i ~N (p y i , σI p ) where p j is a cluster center, σ a given standard deviation and I p the identity matrix of size p. We set one cluster center at p 1 = (1, 1, . . ., 1) and the second cluster center at p -1 = (-1, -1, . . . , -1). In Figures 1A-D respectively, we show four plots of the simulated data reduced to two dimensions thanks to a UMAP [START_REF] Mcinnes | UMAP: Uniform Manifold Approximation and Projection[END_REF], with standard deviations set to 2, 6, 10 and 14 respectively. Naturally, when the standard deviation increases the data becomes less separable.

Model

We apply a simple DNN model, composed of two layers, a first hidden layer with 256 hidden nodes and a classification layer with two nodes, the model is depicted in Figure 1E.

In particular, we minimise a cross-entropy error with an Adam optimiser and a batch size of 16. The tunable parameters are the weight-decay, drop out and learning rate. As an extra regularisation we use batch normalisation.

Application to Histopathology Data

Data Generation and Annotation

The data set used was generated at the Curie Institute and consists of annotated H&E stained histology needle core biopsy sections at 40× magnification sampled from a patient suffering from TNBC. In this paper, we evaluate the prediction of the response to treatment based solely on a biopsy sectioned prior to the treatment. As discussed in the introduction, not all patients respond to NACT, and we are therefore aiming at predicting the response to NACT based on the biopsy. In particular, each section was quality checked by expert histopathologists.

For each patient, we also collect WSI after surgery, allowing an expert pathologist to establish the residual cancer burden, as a proxy for treatment success. Out of the 336 samples that populate our data set, 167 were annotated as RCB-0, 27 as RCB-I, 113 as RCB-II and 29 as RCB-III. This data set is twice as large as the data set used in our previous study (Naylor et al., 2019). Similarly to this study, we refine the number of classes in order to avoid the problem of under-represented class. We investigate two prediction settings: 1) pCR (no residuum) vs. RCB (some residuum) and 2) pCR-RCB-I vs. RCB-II-III, which is clinically more relevant, as it is informative of a patient's prognosis.

Data Encoding

As each biopsy section is relatively big, we wish to reduce the computational burden of feeding the entire biopsy to our algorithms. Instead, given a magnification factor, we divide each biopsy into tiles of equal sizes, 224 × 224 and project this tile into a lower dimensional space. We use a pre-trained DNN on ImageNet [START_REF] Jia | ImageNet: A Large-Scale Hierarchical Image Database[END_REF]) such as ResNet [START_REF] He | Deep Residual Learning for Image Recognition[END_REF] which produces a encoding of size 2048. This process is illustrated in Figure 2 where each biopsy section is converted into an encoded matrix of size n i × P where P is the size of the resulting encoding and n i the number of tiles tissue extracted from tissue i, i ∈ N.

In Table 1 we show the average number of tiles, n i and variance at different magnification factors: highest resolution i.e. no down-sampling (2 0 = 1), down-sampling by a factor 21 and by a factor 2 2 = 4.

The size of the data remains relatively large even after this reduction. We further reduce the size of the tile encoding with a PCA [START_REF] Jolliffe | Principal Component Analysis[END_REF], and project each tile encoding into a space approximately 10× smaller. By keeping 256 components, we keep 93.2% at 2 0 , 94.0% at 2 1 and 94.3% at a magnification factor of 2 2 of the explained variance.

Mathematical Framework

The data set will be denoted by D (X i , Y i ) i∈ [[1,N]] , where every item indexed by i in D is a joint variable (X i , Y i ), N is the size of the data set, X i is the input sample and Y i the corresponding label. As described in the previous section, each tissue is TABLE 2 | Model architecture and pooling strategy, as described in the text. CHOWDER has been published previously [START_REF] Courtiol | Classification and Disease Localization in Histopathology Using Only Global Labels: A Weakly-Supervised Approach[END_REF].

represented by a bag of tiles of variable sizes, in particular ∀i ∈ [(1, N)], X i ∈ R ni×P and Y i ∈ (0, 1) for task (1) or (2). This is simply a multiple instance learning framework, and such a framework has already been implemented for histopathological data [START_REF] Xu | Contextconstrained Multiple Instance Learning for Histopathology Image Segmentation[END_REF][START_REF] Xu | Parallel Multiple Instance Learning for Extremely Large Histopathology Image Analysis[END_REF][START_REF] Courtiol | Classification and Disease Localization in Histopathology Using Only Global Labels: A Weakly-Supervised Approach[END_REF][START_REF] Couture | Multiple Instance Learning for Heterogeneous Images: Training a CNN for Histopathology[END_REF]. We simplify this framework by setting ∀i ∈ i ∈ [(1; N)], n i = n MF 1 which is set accordingly to the chosen magnification factor. For a given sample i, if n i > n MF we down-sample X i , otherwise we up-sample X i to the correct size. We evaluate our models by using the Area Under the Curve of the Receiver Operating Characteristic for measuring the performances in our two binary settings.

Neural Network Architectures

Today, DNN models for WSI classification usually consist in three steps: starting from encodings that are usually provided by pretrained networks, a reduction layer might be applied, followed by an aggregation step that computes a slide level representation from the tile level representations and a final module that maps the slide level representation to the output variable.

In Figure 3, we show a basic example for such an architecture along these lines, with the three algorithmic blocks highlighted in gray. At the tile-level computation, we use 1D convolutions to transform the input encodings X i into a more compact representation X i ′. The tile representations X i ′ are then summarized by a pooling layer, providing us with the biopsy section profile Z i . For this, we can use a standard pooling layer such as an average pooling to quantify the abundance of specific tile patterns, or more complex, attention-based pooling, such as WELDON [START_REF] Durand | WELDON: Weakly Supervised Learning of Deep Convolutional Neural Networks[END_REF]. Finally, from Z i , the slide variable is predicted.

In this article, we test several encoding and agglomeration strategies which are explained in the next sections.

Encoding Projection

In Figure 3, the baseline architecture is illustrated (OneOne), with a 1D convolution for the tile level encoding and one fully connected layer at the slide level. In Figure 4, we test a deeper architecture for the encoding projections, consisting in three consecutive 1D convolutions, including bottleneck layers (depicted in orange), according to best practice in deep learning [START_REF] Huang | Densely Connected Convolutional Networks[END_REF].

Furthermore, we also experiment with skip connections by concatenating the first tile representations to the final representation X i ′ and by concatenating Z i prior to the final softmax. We name this structure ThreeTwoSkip and illustrate it in Figure 5.

Pooling Layers

In terms of pooling layers, we experiment with: average pooling shown in Figure 6A, WELDON [START_REF] Durand | WELDON: Weakly Supervised Learning of Deep Convolutional Neural Networks[END_REF] shown in Figure 6B, a modified version of WELDON shown in Figure 6C and the concatenation of the first and the third is named WELDON-C (for context). The DNN that uses WELDON-C will be named CONAN2 .

The WELDON pooling is an attention-based layer which filters tiles based on a 1D convolution score. In particular, it retains the top and lowest R ∈ N* achieving scores as Z i . This architecture has shown excellent results for specific problems where the biological evidence lies in the detection of one type of specific tiles, like cancer regions [START_REF] Courtiol | Classification and Disease Localization in Histopathology Using Only Global Labels: A Weakly-Supervised Approach[END_REF]. The method however suffers from identifiability issue, i.e., the model can not differentiate between two tiles achieving high or low score. In addition, the agglomeration strategy seems less promising in cases where the information resides in the percentage of tiles of a certain type. By providing a context in which a tile was selected, we allow the model to better differentiate between the selected tiles, thus allowing different tiles with different meanings to be selected, this can be particularly efficient when relevant information is based on different tile patterns.

We recap all the tested models in Table 2.

Baseline Approach

In addition to comparing our proposed architectures to CHOWDER, we also compare them to a much simpler approach where we propagate the slide label to the tile level. If a slide is positive, then we assume all extracted tiles from this slide are positive, if a slide is negative, we assume that all tiles are negative. This is the simplest form of MIL. The training set is therefore huge, and we also known that there will be many errors, as many tiles do not contain any useful information regarding treatment response. Nevertheless, such an approach can work if there is a large fraction of informative tiles.

Model Tuning

We perform a random grid search for most parameters and only in suitable ranges. For the learning rate and weight decay we perform a random log sampling for a random scale associated to a random digit. We range from a scale of 10 -6 to 10 -3 for the learning rate and from 10 -4 to 10 -1 for the weight decay. We randomly sample a drop out from a uniform U [0: 0.4] . We randomly sample a bottleneck layer size from the following list (8,32,64) and the size of the larger representations are randomly sampled from (64, 128).

RESULTS

Simulation Results

High Performance Variability in DNN Training

In Figure 7A, we show the average variance of our model with increasing standard deviation σ. In particular, for each σ, we generate 100 simulated dataset with a standard deviation of σ and train 1000 DNN on the same data and with the same HP. As this is simulated data, we evaluate the performance of each training on a large independently simulated test set instead of using the outer CV loop [START_REF] Varoquaux | Cross-validation Failure: Small Sample Sizes lead to Large Error Bars[END_REF]. We found that setting the learning rate to 1.10 -4 , the weight decay to 5.10 -3 and drop out to 0.4 tend to always return reasonable results for our simulation setting.

As the standard deviation σ of the simulated data increases, we expect more overlapping between our two classes and naturally, the classification accuracy decreases. For lower σ, regardless of using early stopping or not, the models reaches perfect scores.

In Figure 7, we observe that the more difficult the problem (larger σ), the lower the accuracy, but also the larger the variance: not only do we predict less well, but also does the performance variation increase, such that by retraining a model with the same hyperparameters is not guaranteed at all to provide a model with similar performance. We also see from Figure 7A that early stopping alleviates this problem and consistently reduces the variance in performance, in particular for higher σ.

Nested Cross-Validation Leads to Under-performing Models

We compare the performance of the proposed validation procedures to NCV with the number splits k f = 5 in Figure 7B with 95% confidence intervals around the estimator. On the x-axis we have the standard deviation σ of the simulated data and on the y-axis we have the difference between the Naive Bayes estimator with the corresponding performance, therefore the lower the curve the better the procedure. In addition to Figure 7B, we give the corresponding accuracy score for each σ and validation procedure explicitly in Supplementary Table S1. As the dataset distribution is known, the Naive Bayes is the best classification rule that can be implemented and can be viewwed as an upper bound of the performance. For each σ, we collect 100 estimators with Algorithm 1, Algorithm 3 and Algorithm 4. Next, we compare several strategies for NECV and NBFCV: NECV-1/NBFCV-1, where we keep the best scoring model of the ensemble, NECV-3/NBFCV-3, where we keep the top 3 models and NECV/NBFCV, where we keep all five models.

We first notice that the NCV curve is lower or equal to any of the NECV/NBFCV curves. The best performing model is NECV-i.e., the average of all selected models from the inner CV. In particular NECV has a higher Accuracy than NCV by at least 2%. NBFCV under-performs by a slight margin compared to NECV and seems to be equivalent to NECV-3 in terms of performance.

We conclude that retraining the model without an outer validation score leads to lower overall performance and early stopping is a very useful regularization technique for small sample size problems.

Prediction of Response to Neoadjuvant Chemotherapy

We next applied the different architectures detailed in Section 2.4 (also summarized in Table 2) to the problem of the prediction of response to neoadjuvant chemotherapy in TNBC. We tested 3 different image resolutions (0, 1, 2), 0 being the highest resolution. In order to get realistic estimations of the performance, while using early stopping, we perform the validation proposed in Section 2.1. In Figures 8A,B we show the average AUC ROC performance on the residual and prognostic prediction tasks, for all methods shown in Table 2 and for all resolution levels.

For the task of predicting the residual cancer, the best performing model would be the CONAN-c model at resolution 2 with an AUC of 0.654 ± 0.049. Others model's performance range in between 0.55 and 0.60 of AUC with higher standard deviations. Models at resolution 0 seem to generally achieve higher scores then those at lower resolutions. Model architecture c seem to be better suited for this task than the others. The Average concatenated to WELDON-C pooling seems to perform slightly better then the rest. The method CHOWDER which gave excellent results on CAMELYON for cancer detection [START_REF] Courtiol | Classification and Disease Localization in Histopathology Using Only Global Labels: A Weakly-Supervised Approach[END_REF] and which has been a state-of-the-art solution in the field underperforms on our dataset for response prediction.

For the task of predicting the patient prognostic, the best performing model would be the Avg-b model at resolution 0 with an AUC of 0.601 ± 0.019. CONAN-c at resolution 2 performs similarly but with a much higher standard deviation. Neither resolution, nor model architecture and pooling layer seem to unanimously be better then the others. However, CHOWDER under-performs compared to the other proposed methods.

DISCUSSION

In this study, we set out to predict the response to neoadjuvant chemotherapy in TNBC from biopsies taken before treatment. A system that would allow to predict this response with high accuracy could help identifying patients with no or little benefit of the treatment and therefore spare them the heavy burden of the therapy.

From a methodological point of view, this is particularly challenging for three reasons: first, we do not know to which extent the relevant information is actually present in the image data. In addition, even if the relevant information is contained in the slide, the complexity of the related patterns is unclear. Second, biopsies only capture a part of the relevant information, as they are only a localized sample of the tumor. Third, as this is a project regarding a specific subtype, the cohort is relatively small, unlike many pan-cancer cohorts used in large Computational Pathology projects [START_REF] Campanella | Clinical-grade Computational Pathology Using Weakly Supervised Deep Learning on Whole Slide Images[END_REF].

In order to solve this problem, we have developed the model CONAN, that combines the power of selecting K tiles (top and bottom), but keeps both the ranking scores and the full tile descriptions to build the slide representation. We have compared this model with a number of different architectures, and achieved an AUC of 0.65.

We also tackled an important problem of model selection with cross-validation, a crucial step in particular for small datasets. We found that the retraining step in classical Nested cross-validation can lead to lower performances for small N, because the training is highly variable, and a network retrained with the optimal set of hyperparameters is not guaranteed at all to be optimal itself. We therefore have proposed a new cross-validation procedure relying on ensembling rather than retraining, and thus allowing to use early stopping as a regularization method.

Nevertheless, we must conclude that the prediction of treatment response is probably one of the hardest problems in Computational Pathology, and that even though we see that there is some degree of predictability, the results still seem far from clinical applicability. Also, another problem that is not addressed in this study is the applicability of trained networks across different centers. Clearly, we need more data to tackle these challenging questions. But it is also likely that even with much more data, AUCs will not reach very high levels by looking at biopsies alone. A promising avenue would therefore be to use other kinds of data in addition to histopathology data.

FIGURE 1 |

 1 FIGURE 1 | Data set simulation with p = 256 and model, (A) σ = 2, (B) σ = 6, (C) σ = 10, (D) σ = 14. (E) Simple two layer DNN model for predicting the class label.

Algorithm 3 .

 3 Model selection, f ecv .

FIGURE 2 |

 2 FIGURE 2 | Encoding a biopsy.

FIGURE 4 |

 4 FIGURE 4 | Tile computation: Three layers.

FIGURE 5 |

 5 FIGURE 5 | ThreeTwoSkip.

FIGURE 6 |

 6 FIGURE 6 | Aggregation layers: (A) Average pooling; (B) WELDON pooling; (C) WELDON-C pooling which is WELDON concatenated with previous tile encoding.

FIGURE 7 |

 7 FIGURE 7 | (A) Repeated training, with and without early stopping: (A) Performance Variance, as measured on an independently simulated test set, (B)Comparison of models trained with NCV, NECV and NBFCV. On the y-axis we have the difference in accuracy between the Naive Bayes model and a given validation approach. As the data distribution is known, the Naive Bayes is a theoretical upper limit of the achievable performance. Curves are shown with 95% confidence intervals. A lower curve implies a better score.

FIGURE 8 |

 8 FIGURE 8 | Average AUC ROC performance with standard deviation on five fold NECV on the task of predicting: (A) residual cancer and (B) patient prognostic in TNBC.

TABLE 1 |

 1 Mean number of tiles.

	Down-sampling Factor	n i : Mean Number of	i n i
		tiles
	2 0	11,186 ± 6,983	3,758 389
	2 1	2,757 ± 1,757	926,409
	2 2	628 ± 403	211,077

FIGURE 3 | Model OneOne: tile encodings X i from a pretrained network are projected to a more compact representations X i ′ and then aggregated to build the slide representation Z i , which is then used for prediction. Same colors indicate identical dimensions.

Frontiers in Signal Processing | www.frontiersin.org June 2022 | Volume 2 | Article 851809

For Magnification Factor.Frontiers in Signal Processing | www.frontiersin.org June

| Volume 2 | Article 851809

Context cOncatenated tile selection NeurAl Network Frontiers in Signal Processing | www.frontiersin.org June 2022 | Volume 2 | Article 851809

ACKNOWLEDGMENTS

The authors would like to thank André Nicolas from the PathEx platform, Service the Pathologie, Institut Curie, Paris, where the images have been acquired.

DATA AVAILABILITY STATEMENT

The datasets presented in this article are not readily available because they are the property of Institut Curie. Requests to access the datasets should be directed to thomas.walter@mines-paristech.fr.

AUTHOR CONTRIBUTIONS

PN was main contributor, developed the methods and code and generated the results. TL helped with the code. FR and TW designed the project. GB, ML, AVC, A-SH, and FR generated the cohort. GB and ML provided interpretation of the biological results. TW supervised the project. PN and TW wrote the paper.

SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frsip.2022.851809/ full#supplementary-material Conflict of Interest: The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Copyright © 2022 Naylor, Lazard, Bataillon, Laé, Vincent-Salomon, Hamy, Reyal and Walter. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.