Emmanuel Moebel

Antonio Martinez

Damien Larivière

Julio Ortiz

Wolfgang Baumeister

Charles Kervrann

3D ConvNet improves macromolecule localization in 3D cellular cryo-electron tomograms

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The last decades of research in cell biology have revealed that cellular processes are performed by groups of interacting macromolecules in a crowded environment. This is in opposition to previous cell models where macromolecules where considered as isolated objects floating randomly in the cytoplasm.

Deciphering the underlying interaction mechanisms is thus of paramount importance to gain a deeper understanding of the cell. To address this issue, cryo-electron tomography (cryo-ET) is a unique imaging technique capable of producing 3D views of large portions of a cell while having enough resolution to localize and identify macromolecules. Cryo-ET allows avoiding the use of markers, such as fluorescent probes used in light microscopy, which could perturb the cell. First, samples are vitrified in order to preserve both the spatial distribution and the structure of macromolecules in the cell during the image acquisition process. This technique enables the study of cells in a close to native state, and has thus the potential to create a molecular atlas from all the detected components (macromolecules, membranes) observed in cryo-tomograms. However, the analysis of such images is challenging due to poor signal-to-noise ratios and imaging artifacts caused by limited-angle tomography. As a consequence, cryo-ET analysis is heavily dependent on computational tools for interpreting the image content.

A well-established method for localizing macromolecules is template matching (TM) [START_REF] Best | Localization of protein complexes by pattern recognition[END_REF],

where a template containing the macromolecule of interest is used to explore a given 3D cryo-tomogram.

While TM is efficient for localizing large macromoelecules such as ribosomes, it is necessary to apply several image post-processing and analysis methods to decrease the false positive (FP) rate (see Fig. 1). These methods include selection of regions of interest (e.g. areas in the cytoplasm), and thresh-olding the TM score (e.g. correlation score) values. Additional difficulties also arise when TM is used to detect and localize several types of macromolecules that are structurally similar or specific macromolecule states, like binding states of ribosomes (e.g. membrane-bound vs cytoplasmic ribosomes). In general, TM is applied several times to detect each subclass of interest. Unfortunately, the score values are not selective enough to allow one to perform a satisfying classification, especially when the number of considered subclasses is high. Therefore the sub-volumes containing the detected macromolecules (also named particles) of interest are manually analyzed or automatically post-processed by sophisticated classification algorithms [START_REF] Förster | Classification of cryo-electron sub-tomograms using constrained correlation[END_REF]. Such complex and time-consuming processing chains are routinely applied to accurately localize macromolecules and to identify the related native structure in the cell. Note that each TM and sub-volume classification round can each take 10 to 30 hours of computation on specialized CPU clusters.

In this paper, we propose an unified deep learning-based framework [START_REF] Lecun | Deep learning[END_REF] to jointly and fastly localize and classify macromolecules in cryo-ET. Deep learning (DL) is a set of machine learning techniques capable to produce state-of-the-art results in various fields (e.g. computer vision [START_REF] Lecun | Convolutional networks and applications in vision[END_REF], language processing [START_REF] Hinton | Deep neural networks for acoustic modeling in speech recognition[END_REF], super-resolution microscopy [START_REF] Ouyang | Deep learning massively accelerates super-resolution localization microscopy[END_REF] and bioinformatics [?]). In particular, convolutional neural networks (CNN) are able to produce impressive results in image analysis, including image classification [START_REF] Krizhevsky | ImageNet classification with deep convolutional neural networks[END_REF],

segmentation [START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF] and object recognition [START_REF] Szegedy | Deep neural networks for object detection[END_REF]. A neural network is generally composed of successive neuron layers, each transforming incoming data and transferring it to the next layer. The neurons can be seen as small processing units capable of performing linear and nonlinear operations. Each neuron is controlled by parameters which are optimized during the learning process. In the case of CNNs, the neurons are applied in a convolutive manner, which allows dealing with the information redundancy of neighboring pixels (neighboring pixels have similar values). Thus, a neuron can be thought of as a filter, and a neuron layer as a filter bank. The role of a layer is to automatically extract features from the data. Applying sequentially the layers enables to progressively compute more abstract features, which results in a hierarchical representation of the data.

The underlying idea is to learn high-level features from low-level features, which allows a computer to understand complex interactions from basic patterns. The first layers typically encode basic features such as image contours/edges and textures, which allows the next layers to gradually capture more complex shapes (e.g. circles, triangles), objects (e.g. eyes, ears), object ensembles (e.g. faces) and object conditions (e.g. face gender). Those powerful data representations are learned automatically from the data, and tend to be more efficient than conventional handcrafted representations, which require human ressources and are time consuming to design.

Deep learning has been recently investigated to learn high-level generic features in cryo-electron microscopy (cryo-EM). In [START_REF] Wang | DeepPicker : A deep learning approach for fully automated particle picking in cryo-EM[END_REF], the authors proposed first a CNN architecture to au-that the proposed DL framework can be used to pick up ribosomes for subtomogram averaging, an additional post-classification step is necessary to get more satisfying reconstruction results. Unfortunately, no quantitative analysis is presented in [START_REF] Tang | EMAN2 : An extensible image processing suite for electron microscopy[END_REF] to assess the localization accuracy of detected particles and the actual resolution of macromolecule structures once subtomogram averaging is performed. Unlike [START_REF] Wang | DeepPicker : A deep learning approach for fully automated particle picking in cryo-EM[END_REF], Chen et al., 2017], we consider a fully 3D CNN architecture in order to more accurately and reliably detect 3D particles in a native crowded cell environment as illustrated in Fig. 2. In addition, our network is also capable to handle multiple object subclasses at the same time. We especially demonstrate that manipulating jointly a higher number of object subclasses/classes is the key approach to improve performance of CNNs in 3D cryo-ET. Moreover, complex and time-consuming post-classification steps are no longer required to produce reliable results contrary to previous approaches [START_REF] Tang | EMAN2 : An extensible image processing suite for electron microscopy[END_REF]. Besides, while training in [START_REF] Chen | Convolutional neural networks for automated annotation of cellular cryo-electron tomograms[END_REF] is faster (10 min per class) than with our method (12 hours), our processing time is at least twice as fast for one class and remains constant when the number of classes increases (see Fig. 3). In our experiments, we compared our 3D CNN architecture to TM in order to emphasize how the ways in which they were designed and their potentials differ. Unlike TM, our CNN framework is able to localize multiple types of macromolecule at once when applied to a given cryo-tomogram. We also show quantitatively for the first time how the CNN framework, when combined to TM guided procedures, can substantially improve the localization sensitivity of the structure of interest on real cryo-ET data.

The remainder of this paper is organized as follows. In Section 2, we present our deep learning framework. In Section 3, we explain our experimental setup focusing on ribosomes and the results. We show that DL outperforms TM on synthetic images and how TM and DL can be combined to improve the localization sensitivity in real data. In Section 4, we discuss the potential of DL in cryo-ET and future work. clustering algorithm is then applied to aggregate voxels into clusters and to determine the position of particles (gravity center of clusters) in the volume (see Fig. 1). The detected particles are further exploited for subtomogram averaging.

Step #1: Multiclass voxelwise classification

Our objective is to provide a classification map for which each voxel in the 3D map is assigned an object class. The convolutional neural network is usually designed to produce a single output label: it exploits global information by progressively down-sampling the image layer by layer, in order to preserve only relevant information and reduce computation. The disadvantage is that by doing so, the network loses the local information needed for voxelwise classification tasks. While the network is able to reliably decide if an object is present or not in an image, it is not able to accurately estimate the (once training is completed). We compare our method to [START_REF] Chen | Convolutional neural networks for automated annotation of cellular cryo-electron tomograms[END_REF] and to template matching. We do not consider post-processing in displayed time values (i.e. without clustering step for our method, and without connected component analysis and post-classification for [START_REF] Chen | Convolutional neural networks for automated annotation of cellular cryo-electron tomograms[END_REF]). We use a Tesla K80 GPU for our method and a 32-core CPU cluster for template matching, while in [START_REF] Chen | Convolutional neural networks for automated annotation of cellular cryo-electron tomograms[END_REF] the authors used a 12-core CPU workstation.

position of the object. As opposed to this conventional approach, fully convolutional networks (FCNN) [START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF] overcome this difficulty by explicitly combining global and local information in order to provide high resolution label maps. FCNNs involve more interactions in the network and adapt conventional classification networks to perform more robust image voxelwise classification. This idea was exploited in [START_REF] Ronneberger | U-Net: Convolutional networks for biomedical image segmentation[END_REF] and then adapted in [START_REF] Milletari | V-Net : fully convolutional neural networks for volumetric medical image segmentation[END_REF] for 3D image analysis. As described in [START_REF] Ronneberger | U-Net: Convolutional networks for biomedical image segmentation[END_REF], Milletari et al., 2016], our architecture consists of a down-sampling path needed to generate global information and a up-sampling path used to generate high-resolution outputs, i.e. local information (see Fig. 4). Down-sampling is performed with max-pooling layers (factor 2) and up-sampling with up-convolutions [START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF] (sometimes called "backward convolution"), which is basically a trained and non-linear up-sampling operation.

Combining global and local information is performed by concatenating features at different spatial resolutions. The features are then processed with the convolutional layers of the up-sampling path.

Unlike [START_REF] Milletari | V-Net : fully convolutional neural networks for volumetric medical image segmentation[END_REF], our architecture is not so "deep" since we found that using more than two down-sampling stages does not increase the classification results. Also, we used only 3 × 3 × 3 filter sizes as in [START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF]. The rationale behind this choice is that two consecutive 3×3×3 filters mimic a larger 5×5×5 filter but with fewer parameters. Training is then faster and easier and requires less memory. An important concept in neural architectures is the receptive field of deepest neurons layers. It determines the size of the spatial context to be used to make decisions. Considering a large spatial context is essential to handle an object class involving interactions with the environment, for instance interactions with the cell membrane. It is established that adding convolutional layers after down-sampling operations is appropriate to enlarge the spatial context [START_REF] Milletari | V-Net : fully convolutional neural networks for volumetric medical image segmentation[END_REF].

Accordingly, we added two supplementary convolutional layers in the lowest stage of our architecture.

To complete the description, we use rectified linear units (ReLU) [START_REF] Krizhevsky | ImageNet classification with deep convolutional neural networks[END_REF] as activation function for every layer except the last one which uses a soft-max function. While ReLU is a popular choice to tackle non-linearities in the network, the soft-max function is mandatory in order to interpret the network outputs as probabilities for each class.

In summary, our proposed CNN architecture is capable of robustly classifying the cryo-ET tomogram into N subclasses/classes with a high accuracy. Given the voxelwise classification map, the next step consists in estimating the position of each individual object, as described in the next section.

Step #2: Clustering for macromolecule localization

Given the multiclass voxelwise classification map and classification errors, our objective is determine the position of each particle corresponding to a state of a given macromolecule or several types of macromolecules. The voxel labels should be ideally spatially well clustered into well distinct 3D connected components, each cluster corresponding to unique object/particle. Because of noise, non- stationarities in the background, and artifacts in the tomogram, the CNN method generates isolated labels or very small groups of voxels, and groups that contain different label types. Post-processing is then necessary to assign an object class to a given position.

To address this issue, we apply a basic clustering approach. The clusters are built by aggregating neighboring voxels into objects and form 3D connected components. The smallest clusters are considered as false positives and are discarded. The cluster centroid is further used to estimate the object position. As the centroid is computed by uniformly averaging the coordinates of cluster voxels, we are able to numerically produce positions with sub-voxel precision. As several voxel labels can be spatially grouped in a given cluster, the most frequent label or subclass/class is assigned to the detected particle. To address the computational issues, we used the popular mean-shift clustering algorithm [START_REF] Comaniciu | Mean Shift : a robust approach toward feature space analysis[END_REF]. The main advantage of mean-shift is that it is controlled by only one parameter, commonly called the bandwidth, which is directly related to the average object size.

The K-means algorithm was not considered further since the number of clusters must be provided as input parameter by the user.

Training

In training step, the CNN is learned from pairs of tomograms and their corresponding voxelwise classification. In other words, the CNN needs every tomogram voxel annotated as member of the class of interest or as background. While voxelwise classification examples are naturally available for synthetic data, it is often not the case for experimental data. In our case, the experts accurately localized the macromolecules of interest in several training tomograms. Actually, voxelwise classification cannot be performed manually for two reasons: i/ it is time consuming to label each voxel in hundreds of objects in 3D; ii/ the data is so noisy that the object borders are barely visible. To address this issue, we propose an original computational approach based on subtomogram averaging [START_REF] Förster | Structure determination In Situ by averaging of tomograms[END_REF] to label voxels, only from the expert annotations corresponding to the spatial coordinates of macromolecules (see Fig. 5). Subtomogram averaging is a registration algorithm designed to obtain higher resolution structures by averaging thousands of aligned subvolumes containing the same structural unit. The labeled coordinates serve here as inputs to a subtomogram averaging procedure. The subtomograms around the annotated positions are extracted, aligned and finally averaged. The alignment procedure outputs the object orientations, whereas the averaging process provides a clean and missing wedge free density of the macromolecule. From this density, it is possible to create a binary mask of the macromolecule by thresholding the averaged subtomogram. Furthermore, the resulting 3D mask is pasted into an empty volume at each labeled position with the estimated 3D orientation. The resulting volume with well delineated macromolecules is then used as a target to train the parameters of the CNN architecture. It is worth noting that annotating the macromolecule with this semi-automatic approach saves time but may introduce "label noise" in the training. Indeed, we use an average shape to label the macromolecule, and we neglect structural macromolecule variability mainly localized in the object borders. Nonetheless, it has been shown that CNNs have a natural robustness to reasonable amount of "label noise" [Rolnick et al., 2017], which is also confirmed in our experiments.

Due to memory limitations, it is not feasible to load the whole tomogram set with the corresponding targets during training. Therefore, we randomly draw smaller 3D patches around macromolecules at each training iteration. The patch size should be large enough to capture sufficient context information; the macromolecule radius being 10 voxels, we choose a patch size of 56 × 56 × 56× voxels. It is also common to use "data augmentation" when training a CNN; it allows to increase the training set artificially by applying geometric transform to the training images. In our approach, we implement "data augmentation" by applying a 180 • rotation w.r.t. the microscope tilt-axis to each training example. Nevertheless, we do not use typical mirror operations or geometric deformations because the structure of expected objects is the principal clue in the detection problem. Also, we do not use random rotations because of the well-determined orientation of missing wedge artifacts, which is preserved when applying 180 • rotations w.r.t. the tilt-axis. In our experiments, the CNN has been computationally trained for 6000 iterations with the ADAM algorithm, chosen for its good convergence rate [START_REF] Kingma | Adam: a method for stochastic optimization[END_REF], using 0.0001 as learning rate, 0.9 as exponential decay rate for the first moment estimate and 0.999 for the second moment estimate. We use categorical cross-entropy as a loss function (see Fig. 6). The training has been performed on a Nvidia Tesla K80 GPU and took 12 hours of computation, which is reasonable knowing that for other tasks, CNN training can last several days [Krizhevsky et al., 2012] [Simonyan and[START_REF] Simonyan | [END_REF]]. In the next section, we will present the training datasets and the results of our CNN approach applied to both synthetic and real tomograms.

Results

The method has been evaluated and compared to TM on a synthetic and real datasets described below.

Description of data

Dataset #1: synthetic data

The data has been generated with the AV3 toolbox [START_REF] Förster | Structure determination In Situ by averaging of tomograms[END_REF], using atomic densities from the PDB databank [rcsb.org]. Simulation parameters include a voxel size of 13.68 Å(Angströms), a defocus of -6µm, and several values of signal-to-noise (SNR) from 0.05 to 0.10, and tilt ranges from ±50 • to ±70 • , with a tilt-increment of 2 • . In this experiment, nine classes of prokaryotic macromolecules have been chosen to depict varying amounts of inter-class similarity (see Fig. 7). The GroEL and ribosome classes are significantly different in size and shape, whereas different functional states of the proteasome and GroEL are structurally similar. In addition to these well-known macromolecules, we have introduced basic objects (ellipses, spheres, discs) in order to mimic not well-defined structures found in a cell, like membrane components, small molecules, and gold particles.

The training set consists of 510 macromolecules for each class, the validation set is composed of 240 macromolecules for each class, and the test set is made of 105 macromolecules for each class. The macromolecules and basic objects have been placed at random positions and orientations in the 3D volumes in order to simulate the crowded environment of a cell.

Dataset #2: experimental data

The second dataset is composed of 63 tomograms of Chlamydomonas Reinhardtii cells (see [START_REF] Pfeffer | Dissecting the molecular organization of the translocon-associated protein complex[END_REF] for details about data acquisition) and has been annotated for membrane-bound 80S ribosomes (denoted mb-ribo in the following) positions by experts. To get these annotations, the experts first used TM with a template generated from the dataset, using manually selected ribosomes and subtomogram averaging. Then they refined the TM results by applying subtomogram classification (CPCA [START_REF] Förster | Classification of cryo-electron sub-tomograms using constrained correlation[END_REF]) and performing careful visual inspection (see Fig. 1). To reduce computational cost, the tomograms were under-sampled, resulting into a tomogram size of 928 × 928 × 464 voxels and a voxel size of 13.68 Å. Tilt range is ±60 • with an increment of 2 • . In the end, the dataset has been annotated with 9487 mb-ribos. The subtomogram average computed from the total set of mb-ribos is the best model we have for 80S ribosomes in Chlamydomonas Reinhardtii cells, and is therefore used as a reference for subtomogram alignment. As this dataset was originally annotated and designed to study the (mb-ribo) class, we used computational tools to get examples corresponding to the cytoplasmic ribosome (denoted ct-ribos) class and the membrane class.

First we added the membrane class by employing an algorithm dedicated to membrane segmentation [START_REF] Martinez-Sanchez | Robust membrane detection based on tensor voting for electron tomography[END_REF]. Meanwhile, we got (ct-ribo) examples by applying TM and selecting the most isolated candidates, located at a distance higher than 273.6 Å (i.e. the ribosome diameter) to membrane components. The motivation behind adding new classes to the available annotations was twofold: on the one hand, our motivation was to demonstrate the multiclass ability of our method on real data; on the other hand, we noticed that the multiclass approach tends to improve the discriminating power of the network and the capacity to reliably detect mb-ribos in real tomograms.

When trained by using only the mb-ribo examples, the network actually detects unwanted cytoplasmic ribosomes. By considering the ct-ribo class in the training, we encourage the network to better discriminate both ribosome subclasses (corresponding to binding states). Note that these annotations of membrane components and cytoplasmisc ribosomes have been obtained without the supervision of an expert. Therefore more errors are expected for these two classes when compared to the mb-ribo examples reliably annotated by the experts in our protocol.

Evaluation

The metric used to assess localization performance is the F1-score, as commonly used in detection problems. The F1-score can be interpreted as a weighted average of two well-known metrics depending on the number of true positives (TP):

• Recall (also known as sensitivity):

R = Number of TP Number of particles in the tomogram

,

• Precision (also known as positive predictive value):

P = Number of TP Number of localized particles . (2)
The F 1 score defined as follows

F 1 = 2 RP R + P ,
allows one to evaluate performance by considering a single value.

In what follows, a detected particle is considered as a TP if is closer than 136.8 Å (i.e. ribosome radius) to a ground truth object.

First, we quantitatively evaluated the performance of the multiclass localization method on dataset #1. Synthetic data are very helpful to objectively study the influence of SNR and tilt-ranges on results.

To perform multiclass localization with TM (which is mono-class), TM has been applied once for each class with the corresponding templates. The templates have been computed from the atomic densities used to generate the dataset #1. A threshold applied to the TM score map is used to select the N best candidates for each class. The threshold is chosen so that the overlap with the ground truth is maximized. As recommended in [START_REF] Best | Localization of protein complexes by pattern recognition[END_REF], if multiple templates compete for the same position, the class with the highest score is selected.

Dataset #2 was used to evaluate the performance of our method in real conditions. In our protocol, only the mb-ribo annotations have been provided by experts, and are therefore considered as the ground truth. Accordingly, we only quantify results for the mb-ribo class, while providing visual results for the ct-ribo and membrane classes. Unlike [START_REF] Tang | EMAN2 : An extensible image processing suite for electron microscopy[END_REF], we also analyzed the score distributions and subtomogram averages. The resolution of obtained subtomogram averages is estimated with the commonly-used gold standard Fourier shell correlation (FSC) combined to the "0.143" threshold criterion.

Result analysis

Dataset #1: synthetic data

It turns out that DL outperforms TM for each class, in terms of F 1 -score (see Fig. 7 A). An interesting result is that the scores achieved by DL are virtually perfect. In comparison TM scores are class dependent, ranging from 0.07 to 0.98. TM achieves good scores for the biggest macromolecules (ribosome 70S and FAS), but the performance is lower if the size of the macromolecule is small. Indeed, the confusion matrices (see Fig. 7 B) reveal that small targets like proteasomes and GroEL are often confused with background or decoy objects. In addition, TM has some difficulty to tackle inter-class similarity, especially for the three functional states of the proteasome (double bpa, single bpa and without bpa).

We evaluated the robustness of TM and DL methods by varying the signal-to-noise ratio (SNR)

and tilt-range. Figure 8 plots the average F 1 -score over the nine classes. As expected, we observe a performance drop for decreasing values of SNR and tilt-range. However, DL is remarkably stable and produces nearly constant scores. TM on the other hand loses 6% of F 1 -score when the SNR decreases from 0.15 to 0.05, and 7% when the tilt-range decreases from ±70 • to ±50 • . In all tested situations, it turns out that DL is more robust to noise and missing-wedge than TM. In summary, the results on dataset #1 prove that DL is capable to provide better results than TM as soon as DL exploits ideal ground truth during training, that is with no "label noise". These results demonstrate that DL outperforms TM on synthetic noisy and MW corrupted data.

Notice that our DL approach is able to implicitly tackle the missing wedge (MW), while TM usually considers the so-called constrained cross correlation to handle the MW information. Actually, DL is capable to manage non-linear object deformations due to MW during training, without additional prior imposed by the experts.

Dataset #2: experimental data

Comparison of score values of TM and DL In the first part of experiments, we have carefully examined the scores produced by TM and DL methods. In Fig. 9, it is clear that the TM score map is much noisier than the maps generated by DL. Actually, the responses of TM based on the constrained cross correlation are very high at ribosome locations and at undesirable locations corresponding to highly contrasted structures with similar sizes (for instance see cell membrane in Fig. 9). TM tends to generate a lot of false positives in the cell. Consequently the experts need to apply post-processing techniques to select relevant information in order to exploit the TM results. Unlike TM, DL provides clean score maps and only depict high values in well-localized blobs. These results suggest that DL is more capable to properly learn the structure and geometry of complex macromolecules.

In order to further support this idea, we examined the distribution of local maxima values in each score map (see Fig. 9). A sharp mode (depicted in red in Fig. 9) can be regarded as an indicator to assess discrimination quality. As shown in Fig. 9, the mode for TM is weak, while for the DL ct-ribo class the mode is not very sharp either. However, for the DL mb-ribo class, the mode is more significant and sharp, suggesting a score less prone to ambiguity. This observation confirms the idea that the TM score is actually not very discriminating in general. More specifically, it is not a surprise if the score of the ct-ribo class is not as discriminating as the score of the mb-ribo class. This is probably related to the annotation quality used for learning, much more higher in the case of mb-ribo class. In our experiments, the mb-ribos annotations have been carefully performed by an expert.

In summary, DL produces more sharper scores than TM (constrained cross-correlation), when the training data is carefully labeled (we have a better performance for mb-ribos than for ct-ribos). the Euclidean distances between the ribosomes and the nearest membrane component (see Fig. 10 B).

Evaluation of voxelwise multiclass classification

The distance histograms confirm the visual analysis: for mb-ribos, the distance histogram has a sharp mode at 136.8 Å, which corresponds to the ribosome radius. Therefore a large majority of the voxels classified as mb-ribos are very close or linked to the membrane. For ct-ribos, we notice no sharp mode and the distance histogram looks Gaussian with a large variance. This confirms the heterogeneity of ribosomes freely floating in the cytoplasm. In conclusion, our method allows accurate multiclass object detection in cryo-ET.

Evaluation of overlap with annotations The next step of our evaluation process consists in measuring the overlap between the expert annotations and the mb-ribos found by DL as explained earlier.

We compare the results to the TM outputs, that is before applying sophisticated post-classification methods. In Fig. 11, we plotted the Recall, Precision and F 1 score w.r.t. the DL and TM algorithm parameters. We focused on the thresholds used to detect objects (TM: threshold on score values; DL: object size threshold) (see Sec. 2.1.2). We obtained a F 1 score of 0.86 for DL and a F 1 score of 0.50 for TM. These numbers illustrates the ability of our DL approach to learn and bypass the expert processing chain. Moreover, the computation time of our DL approach is very small when compared to the TM algorithm as given in Figure 3. Now that it has been established that DL has a better overlap with the annotations than TM, in the remaining we focus our analysis DL detections.

We have also examined the complementarity between the two sets of mb-ribo macromolecules detected by the the experts (guided by TM) and the DL method. In what follows, we respectively denote S E and S DL the sets obtained by the experts and the DL method. While the overlap S E ∩ S DL between both sets is substantial (1516 particles), there is also a significant amount of particles belonging

to S E \ S DL (220 particles), i.e. the particles annotated by the expert but overseen by DL, and to S DL \ S E (356 particles), i.e. particles found by DL but overseen by the expert. We can benefit from the two complementary object position estimations to improve overall validation rates. Actually, the union S E ∪ S DL of the two sets enables to increase the list of potential mb-ribo macromolecules, for which a confidence level can be assigned to each member depending on whether it belongs to S E ∩S DL , S DL \ S E or S DL \ S E . Objects belonging to S E ∩ S DL , i.e. found by both methods, are very likely to be true positives. Meanwhile the detected objects belonging to S E \ S DL and S DL \ S E can be labeled as "suspicious" and need more investigation. These two sets are relatively small and the experts may focus on the detected macromolecules that may correspond to rare conformations observed in the cryotomogram. From our analysis, it is possible to get a high overlapping rate with the expert annotations by using our DL approach, suggesting that DL is able to learn the expert analysis chain.

Analysis of subtomogram averaging results

It is usually recommended in cryo-ET to analyse the 3D structure of macromolecules by using subtomogram averaging. Accordingly, we analysed the detected particles by computing subtomogram averages for each ribosome subfamily (see Fig. 12).

In this way, we compare the subtomogram averages obtained from the DL detections and the expert annotations. The same process has been used to compute all averages, namely fast rotational matching [START_REF] Chen | Fast and accurate reference-free alignment of subtomograms[END_REF] (see Sec. 3.4). The averages of mb-ribo are composed of two densities: the average ribosomal density and the average membrane density. The intensity level of the membrane density varies with the proportion of mb-and ct-ribos in the average; the higher the intensity is, the more the number of mb-ribos is high.

As expected, the average membrane density computed from the set S DL of mb-ribo particles has a high intensity, while the average membrane density computed from the set of ct-ribo particles detected by our DL approach, is non-existent. Consequently, our method makes little to no confusion between mb-and ct-ribos. We notice that the average ribosomal density is similar in both ribosome subfamilies, the difference lying in the neighborhood of the link with the membrane. Therefore, DL is able to efficiently represent the geometric structure of the macromolecule, and at the same time to capture the local context and interactions of the macromolecule with the environment.

Next, we computed the Fourier shell correlation (FSC) for obtained subtomogram averages (see Fig. 12). According to Fig. 12, the best resolution has been obtained with the expert average: 24 Å versus 24.7 Å and 32.7 Å for the DL subclasses mb-ribo and ct-ribo respectively. Our method therefore allows to achieve a resolution comparable to an expert. As for the ct-ribo average the resolution is lower, most likely because the ct-ribo annotations used for learning are of lower quality (see Sec. 3.1.2).

Even though resolutions for the expert and DL mb-ribo averages are very close, the poorer resolution values of the DL averages can be caused by two main factors:

H 1 : The averages contain potential false positives, suggesting that our method probably made a few mistakes.

H 2 : The averages contain particles with a low quality, suggesting that our method found supplementary particles, missed or discarded during the annotation process.

In what follows, we performed further investigations to check the two hypotheses H 1 and H 2 .

First, we decided to align and average the mb-ribo particles of the set S DL \ S E , i.e. found by DL but absent in the expert annotations (see Fig. 13). If no clear signature of ribosome density appears in the average, hypothesis H 1 is valid. On the contrary, if we observe a ribosome patterns in the average, we can conclude that our DL method found additional ribosome particles potentially discarded by the experts (hypothesis H 2 ). Nevertheless, note that the two hypotheses are not mutually exclusive.

In Fig. 13, we display the resulting subtomogram average denoted A DL computed from 356 detected particles. Since it is not guaranteed that all the particles involved in the average are actual mb-ribos, we evaluated the difference between A DL and the average A ref computed from 356 true mbribos (expert annotations) randomly picked from the the set S E ∩ S DL . Also, in order to check if A DL is not biased by the reference template used for subtomogram alignment (see Sec. 3.1.2) as described in [Henderson, 2013], we computed another average denoted A • DL from 356 subtomograms picked from random positions. We compared A DL , A ref DL , and A • DL visually and by estimating the underlying resolution (see Fig. 13). We notice that A DL has a lower resolution (36.4 Å) than A ref DL (33.5 Å). It means that A DL probably contains false positives and/or very noisy instances. Nonetheless, A DL has a higher resolution than A • DL (48.6 Å), suggesting that the reference template bias is not significant. Moreover, a ribosome pattern visually appears in A DL . This suggests that our DL approach has actually found mb-ribos that have been missed during the annotation process.

To be fair, we applied a similar comparison to S E \ S DL , i.e. the set of mb-ribos annotated by This set of results emphasizes that TM and DL can be combined to better investigate cryotomograms. The set of common objects found by the two methods enables to focus on detections exclusively found by the experts or by DL. In addition, we show that, while the FSC curves obtained with DL are below the curve obtained with expert annotations (see Fig. 12), it is risky to decide that the detected macromolecules are not ribosomes (see Fig. 13). The estimated resolutions are lower mainly because the particles involved in the subtomogram averaging are corrupted by noise and other sources of signal degradation. The supplementary noisy particles need to be further examined by experts since they may be valuable mb-ribo candidates. Finally, it appears that the number of actual mb-ribos is higher than expected: in our test set, we have detected +20.5% of mb-ribos when compared to the S E set. In summary, DL found additional noisy mb-ribos that were missed or discarded during the annotation process.

Implementation details of the DL (3D CNN) software

To implement our 3D CCN method, we used Keras [keras.io], an open-source toolbox written in python and using the Tensorflow framework.

As to template matching and subtomogram averaging, we used the PyTom toolbox [START_REF] Hrabe | PyTom: A python-based toolbox for localization of macromolecules in cryo-electron tomograms and subtomogram analysis[END_REF].

We used the in Pytom implemented fast rotational matching routine [START_REF] Chen | Fast and accurate reference-free alignment of subtomograms[END_REF] for subtomogram alignment. The alignment has been performed with respect to a reference template (see Section

3.1.2).

For 3D visualizations, we used Chimera [START_REF] Pettersen | UCSF Chimera -A visualization system for exploratory research and analysis[END_REF].

Discussion

We proposed a 3D CNN framework for voxelwise multiclass classification and particle localization in 3D cryo-ET. We showed on synthetic data that our DL method has superior localization performance than TM. On real data, our DL method is able to discriminate two binding states of ribosomes and to segment cell membranes. We achieved 86% of overlap with expert annotations. Also, when applying subtomogram averaging to the detections, we obtain a resolution comparable to the expert. The supplementary membrane ribosomes found by DL and missed during the annotation process, can be further inspected since the related set is small. While it is established that TM is widely applied in cryo-ET to successfully detect various macromolecular complexes, as of now it has mainly been used to detect relatively large macromolecules like ribosomes. It is worth noting that DL found ribosomes that where missed by the expert annotation process (assisted by TM), and conversely. Consequently both methods can be combined to investigate more efficiently cryo-tomograms. In particular, more effort can be made to analyze the particles detected by either DL or TM, that is when the algorithms are not in agreement. TM and DL can also be combined in a virtuous manner as follows: experts can perform a first round of detection with the routine tools, then train our method from the resulting detections; in a second round, our DL framework can be applied to check the missed objects and increase the size of databases.

To our knowledge, TM is not usually applied to detect more than two or three classes/subclasses on the same dataset. However, cryo-electron tomograms contain much more information. They offer 3D views of the whole macromolecular environment, albeit hardly discernable from noise and artifacts.

Therefore more powerful pattern recognition and machine learning techniques are needed to analyse contents and extract information. DL has been shown to be able to handle the spectacular amount of 1000 classes [START_REF] Krizhevsky | ImageNet classification with deep convolutional neural networks[END_REF]. Also, it is invariant to diffeomorphisms [Mallat, 2016] and then able to cope with non-rigid, elastic deformations within a class. A major disadvantage of TM is that it can only handle rigid views of an object, which is problematic knowing that many proteins (e.g. proteasomes) have a high structural variability. Thus, provided that DL has been trained with enough representative examples capturing shape variability, the CNN should be able to learn different conformations of the same macromolecule.

years, the Fourmentin-Guilbert Foundation has supported the Max Planck Institute of Biochemistry, headed by Wolfgang Baumeister whose team has been capable of delivering whole cryotomograms of E. coli cells at an unprecedented resolution.

The next big step, still an ongoing challenge, was to recognize macromolecular components within the tomograms. Most of the effort of the scientific community was put on the identification of the ribosomes thanks to a methodology relying on single-particle analysis and template matching and giving impressive results. However, it is likely that such an approach will be limited to "big" complexes like the ribosomes.

Facing this challenge, the Fourmentin-Guilbert Foundation has solicited the research group headed by Charles Kervrann to develop alternative recognition methods having the potential to help the in-situ identification of the thousands of proteins left in the dark. The Serpico Proejct-Team has then developped and compared with well-established methods new approaches based on deep learning and capable of identifying and counting the "gold standard" ribosomes within a tomogram. Theses methods, as an alternative to template matching, should also have the potential to apply on particles smaller and rounder than ribosomes.

Figure 1 :

 1 Figure 1: Comparing workflows: on the left, our DL framework and on the right, the common processing chain involving TM. Our approach is multi-class, whereas the TM processing chain needs to be applied once for each class.

Figure 2 :

 2 Figure 2: Chlamydomonas cell. (A) Tomogram slice; blue arrows indicate mb-ribos; yellow arrows indicate ct-ribos. (B) Corresponding voxelwise classification obtained by our 3D CNN, performed for 3 classes: mb-ribos (blue), ct-ribos (yellow) and membrane (gray).

  multiple objects with a CNN architecture Given a training set of object classes, we propose a supervised 3D CNN-based method to classify the 3D tomogram voxels into several types of macromolecules or states of a given macromolecule. A

Figure 3 :

 3 Figure 3: Time needed to process a tomogram of size 928 × 928 × 464 voxels, w.r.t. the number of classes

Figure 4 :

 4 Figure 4: CNN architecture. In green convolutional layers labeled with (#filters × (filter size)). In the last layer, Ncl stands for the number of classes.

Figure 5 :

 5 Figure 5: CNN training: this figure illustrates how to obtain voxelwise classification examples for training, using only position annotations.

Figure 6 :

 6 Figure 6: Evolution of loss and accuracy during training on dataset #2. These quantities are computed for the train set, as well as for the validation set, in order to estimate the generalization capabilities of our network. The curves for both sets should overlap, else it indicates overfitting (the network memorizes train samples instead of learning discriminating features).

Dataset # 2

 2 has been arbitrarily split into training, validation and test sets. Training and validation sets have been sampled from 55 tomograms and consist of 5971 mb-ribos for training and 1493 mb-ribos for validation. The test set is composed of 8 tomograms annotated with 1736 mb-ribos.

  Figure 7: Dataset #1: comparing DL and TM performances per class. For this result we used SNR=0.1 and a tilt range of ±60 • . (A) displays the achieved F1-scores and (B) are obtained confusion matrices, illustrating the miss-classifications of both methods.

Figure 8 :

 8 Figure 8: Dataset #1: average F1-score for each method, for varying SNR (A) and tilt-ranges (B). For (A), we consider a ±60 • tilt-range and SNR values 0.15, 0.10 and 0.05. For (B), we consider a SNR of 0.10 and tilt-ranges ±70 • , ±60 • and ±50 • . The images below the curves illustrate the effects of the varying parameters on a synthetic data sample (in image domain for (A), in spectral domain for (B)).

Figure 9 :

 9 Figure 9: Dataset #2: comparing the score maps obtained from TM and our DL approach. On the score-maps bottom, zoomed-in windows and histograms of local maxima values.

Figure 10 :

 10 Figure 10: (A) 3D voxelwise classification of experimental tomograms, as obtained by our CNN. The classification displays cell membrane (in gray), membrane-bound ribosomes (blue), and cytoplasmic ribosomes (yellow). (B) Distance to membrane histograms of detected ribosomes, on the left for mb-ribos and on the right for ct-ribos.

Figure 11 :

 11 Figure 11: Overlap with expert annotation w.r.t. method threshold parameter: on the left for TM, on the right for our DL approach.

  the expert but missed by DL. As before, we obtain A E from the 220 objects belonging to S E \ S DL , A ref DL by randomly sampling 220 mb-ribos from S E ∩ S DL , and A • DL from 220 random positions. The obtained resolutions are very similar to what is achieved with DL. Here again, A E has a lower resolution (37.6 Å ) than the reference A ref E (34.2 Å ). It is noteworthy that in both cases, the mb-ribos

Figure 12 :

 12 Figure 12: On top, the subtomograms obtained from, the expert annotations, the DL detections for mbribos and ct-ribos, respectively. All averages have been obtained with the same alignment procedure and parameters. For visualization purpose, the averages have been low-pass filtered at 40A resolution. At the bottom, the corresponding gold-standard FSC curves with estimated resolutions.

Figure 13 :

 13 Figure 13: Top-middle: diagram representing the overlap between mb-ribo sets SDL and SE. The emanating arrows represent from which sub-set the displayed subtomogram averages originate. Bottom: FSC curves for each subtomogram average.

Acknowledgment

This work was jointly supported by Fourmentin-Guilbert Scientific Foundation, and Région Bretagne (Brittany Council). Experiments on real data (courtesy of MPI Biochemistry, Martinsried, Germany) were performed on the Inria Rennes computing grid facilities partly funded by France-BioImaging infrastructure (French National Research Agency -ANR-10-INBS-04-07, "Investments for the future").

Note

Through the LifeExplorer project, the Fourmentin-Guilbert Foundation has been a pioneer in approaching the structural modeling and visualization of entire cells. It is expected from the creation of interactive 3D avatars of cellular environments, bridging from the level of atoms to the level of cells, that the rules governing the spatiotemporal organization of the cytoplasm could be revealed. Such an approach requires to make an inventory and a cartography of all the components constituting a single cell.

The technique of choice for such a mapping is cryoelectron microscopy applied on frozen but intact cells. For