Filippo Lunghini

Gilles Marcou
email: g.marcou@unistra.fr

Philippe Azam

Marie-Hélène Enrici

Erik Van Miert

Alexandre Varnek
email: varnek@unistra.fr

Visualization and analysis of the REACH-chemical space with Generative Topographic Mapping

Keywords: REACH chemical space, Generative Topographic Mapping (GTM), environmental fate, ecotoxicology, visualization

In the framework of REACH (Registration Evaluation Authorization and restriction of Chemicals) regulation, industries have generated and reported a huge amount of (eco)toxicological data on substance produced or imported in Europe. The registration procedure initiated the creation of a large REACH database of well defined (eco)toxicological properties. Here, the data distribution in the REACH chemical space was analyzed with the help of the Generative Topographic Mapping (GTM) approach. GTM generates 2-dimensional maps on which each compound is represented as a data point. The 3 rd dimension can be used in order to display a distribution of the given (eco)toxicological property, which can further be used for property assessment of new compounds projected on the map.

We report the "Universal REACH map" which accommodates 11 endpoints, covering environmental fate and (eco)toxicological properties. This map demonstrates acceptable predictive performance: in crossvalidation, balanced accuracy ranges from 0.60 to 0.78. The 11 endpoints profile has been computed for each REACHregistered substance. Some concerns related to acute aquatic toxicity have been identified, whereas for environmental fate and human health endpoints the amount of compounds predicted as of concern was much smaller. It has been demonstrated that superposition of several class landscapes allows to select the zones in the chemical space populated by compounds with a given (eco)toxicological profile.

Introduction

The European REACH (Registration Evaluation Authorization and restriction of Chemicals) regulation [1] , established in 2007, requires from industry to register all substances imported or manufactured in quantities larger than one ton/year. To this end, the registrants must submit to the European Chemicals Agency (ECHA) a technical dossier that characterizes for a given substance, its physical-chemical, environmental fate and (eco)toxicological properties, called endpoints. In order to decrease a need for expensive experimental tests, the REACH regulation allows to use some alternative methods, including predictive statistical models.

In the past years, several models predicting REACHrelevant endpoints were obtained with the help of various machine-learning methods like multiple linear regression, support vector machine or neural networks [2][3][4][START_REF] Echa | Practical Guide How to Use and Report (Q)SARs[END_REF] . However, some of them suffered from absence of technical documentation complying with the REACH requirements [START_REF] Oecd | Guidance Document on the Validation of (Quantitative) Structure Activity Relationship [(Q)SAR] Models[END_REF] , for instance concerning the model's Applicability Domain (AD) or its validation procedure [3,[START_REF] Golbamaki | [END_REF] . In our recent studies [2,4] , we found that for some endpoints, currently-existing tools have disappointing performances when applied to compounds coming from an industrial context, due to their restricted AD. Some existing tools used for the chemicals ranking according to their environmental and toxicological concern, such as DART (Decision Analysis by Ranking Techniques) [8] , use experimental data as an input. To overcome these drawbacks, here, we propose an integrated profiling methodology based on Generative Topographic Mapping (GTM) [9] . GTM is a probability-based mapping strategy, which can be applied both for large-scale data visualization and property prediction. We chose this method because of the following reasons: (i) GTM allows multi-task learning, since several properties can simultaneously be accounted for; (ii) it produces a graphical output (i.e. a two-dimensional map of the chemical space) which can be used as support to better understand the model's output, in the light of a mechanistic interpretation. Here, a set of curated data for 11 endpoints prepared in our previous studies [2][3][4] ( "Global dataset") has been used to train the GTM model and to delineate the REACH chemical space. The following 11 endpoints were considered: bioconcentration factor, ready biodegradability, environmental persistence in sediment, soil and water media, acute aquatic toxicity to algae, daphnia and fish, rat acute toxicity, androgen and estrogen receptor binding potential. Our goal was to demonstrate a potential of GTM to integrate available data into a multi-task modelling framework and to facilitate identification of the compounds of potential concern. For this purpose, a profile assembling selected 11 properties has been computed for each substance registered under the REACH Regulation (called "REACH-INV").

The manuscript is organized in the following way: first, we describe the GTM model built on the Global dataset, followed its application to profile the REACH-INV compounds. Second, we discuss performances of the GTM approach as multi-task learning method, highlighting its potential as a profiling tool.

Materials and methods

Considered endpoints

Table 1 reports the 11 considered endpoints. A brief description is provided in the following paragraphs.

2.1.1

Bioconcentration Factor (BCF) BCF estimates the tendency for a xenobiotic to concentrate inside living organisms. It is defined as the process of concentration of the chemical from the water phase through non-dietary routes, such as absorption from respiratory surfaces (e.g. lungs/gills) or skin [2] .

Ready biodegradability (RB)

Long term exposure for living organisms to many xenobiotics is dependent on the environmental fate of such chemicals which in turn is highly dependent on their biodegradation. Biodegradability is determined by a multistep procedure. This assessment usually starts with a very stringent first-tier assay, providing a binary classification whether the substance rapidly degrades in the environment, called "ready biodegradability" [4] .

2.1.3

Sediment, Soil and Water Persistence (SedP, SoilP, WatP) Unlike relatively cheap and fast ready biodegradability assay, these higher-tier simulation studies are carried out when the substance's degradation half-life (in a given environmental compartment) value actually needs to be evaluated [10] .

2.1.4

Aquatic acute toxicity to Algae, Fish and Daphnia (AlgaeTox, DaphniaTox and FishTox) Aquatic acute toxicity tests aim to estimate the short-term toxicity [START_REF]Guidance on Information Requirements and Chemical Safety Assessment Chapter R.7b: Endpoint Specific Guidance[END_REF] against three species belonging to different trophic levels, considered to be representative of the aquatic ecosystem. Briefly, the test organisms are exposed to the studied substance via the water media, and the following substance-induced effects are measured: (i) for Algae, growth inhibition, expressed as median effective concentration (EC50) measured at 72 hours; (ii) for Daphnia, immobilization at 48 hours expressed as median effective concentration (EC50); (iii) for Fish, the median lethal concentration at 96 hours (LC50).

2.1.5

Rat acute toxicity (RatTox) Rat acute toxicity estimates the short-term lethality (hazard) to humans following ingestion for which oral administration to rodents is used as a proxy. The REACH regulation requires its assessment even for small tonnages. Consequently, this experimental test is one of the most commonly performed animal tests which partly explains its much higher data availability compared to the other endpoints [3] .

2.1.6

Androgen and Estrogen receptor binding (AR binding and ER binding) An endocrine disrupting chemical is an exogenous substance that alters the functions of the endocrine system to the point of causing adverse effects. Possible ways for a chemical to alter the endocrine system is to bind to androgen or estrogen receptors in an agonist or antagonist way. In the framework of the "Collaborative Estrogen Receptor Activity Prediction Project" (CERAPP) [START_REF] Mansouri | [END_REF] and "Collaborative Modelling Project for Androgen Receptor Activity" (CoMPARA) [13] international workgroups, a large number of compounds were tested for their potency to disrupt the AR/ER signaling pathway chains.

Data collection and curation

Experimental data was collected from multiple publicly available databases and scientific literature [2][3][4] . Among them, the main source was the database of the European Chemical Agency (ECHA) [START_REF]eChemPortal: Global Portal to Information on Chemical Substances[END_REF] , which comprises the REACH-registered substances. Raw data processing and standardization were done with workflows implemented in the Konstanz Information Miner (KNIME) software [START_REF] Berthold | [END_REF] . The PubChem [START_REF]The PubChem Project[END_REF] online service was queried to verify SMILES correctness. Generated SMILES were then standardized with the following rules: removal of salts/solvents, removal of explicit hydrogens, aromatic representation of benzene rings, removal of stereo information and transformation of -nitro and -sulpho containing groups into canonical notation, neutralization. Duplicates were removed based on standardized SMILES matching.

As some endpoints are typically described by continuous values (e.g. acute toxicity) while others by categorical values (e.g. ready biodegradability), all the former properties were discretized into "Concern" (C) or "non-Concern" (nC) binary classes (Table 1). For this purpose, REACH-relevant threshold values were selected: for instance, a substance is defined as persistent in sediment (class C) if its degradation half-life is higher than 120 days.

Whenever it was possible, the thresholds have been selected with the goal to make them as conservative as possible. For acute aquatic toxicity, we introduced a ten-fold increase to the 1 mg/L limit (mentioned by REACH to classify the substance as hazardous for the aquatic environment), raising it to 10 mg/L. For the bioaccumulation factor, we selected the cutoff of BCF > 3 log unit, which corresponds to the bioaccumulation limit as defined by the US EPA (that is more conversative than the one reported by REACH of 3.3 log unit). The only exception was for acute rodent toxicity: in that case, a less conservative value of 300 mg/kg (corresponding to "acute toxicity category 4" under the Classification, Labelling and Packaging -CLP Regulation) has been used. To be conversative, the cutoff of >2000 mg/kg (which corresponds to "not classified" under the CLP) should have been selected to discriminate C vs. nC classes. However, this lead to an unbalanced dataset (with the majority of compounds being classified as of concern) leading to less contrasted map that we considered less useful maps in the current situation.

It is important to emphasize that these thresholds have been selected with a certain degree of arbitrariness and their relevance vary according to the field of application (for instance when considering a different regulatory framework).

For the remaining endpoints, the label C was assigned when the experimental values had an assignment of "concern" (e.g. binding to the AR/ER receptors or not readily biodegradable), as opposed to the nC label.

The Global dataset results from the merging of all the available data on the abovementioned 11 endpoints: it counts 29433 experimentally measured data points for 17762 unique compounds (as one compound can have more than one associated experimental value). This dataset has been assembled and curated in our previous studies [2][3][4] . Table 1 reports the sizes of the endpoint subsets. The REACH-INV set comprises the entire inventory of the substances registered under REACH, that have been extracted from the European Chemicals Agency database [START_REF] Echa | ECHA Website[END_REF] , and chemical structures were curated using the same procedure as for the Global set. As this database has been queried in our previous works [2][3][4] , there is a certain degree of overlap between the Global dataset and the REACH-INV set, as it is shown in section 3.

Modelling workflow

The workflow is depicted by Figure 1. Its main steps are described in this section.

Molecular descriptors

The Global dataset, is encoded (Figure 1, step 1) by several types of ISIDA Property-Label Molecular descriptors [START_REF] Ruggiu | [END_REF] . These descriptors work as substructures (fragment) counts of a molecule -for example, D1 = number of "C=O" groups, D2 =number of "C-N-C" fragments, etc. The molecule can be fragmented using two main fragmentation patterns: sequences or atom centered fragment. Moreover, in both cases the size of the fragment (length or radius, respectively) can be varied. Each unique fragmentation scheme is referred as descriptor space. Several tens of descriptor spaces were generated. This list of different fragmentation patterns was used as a starting pool for searching the most appropriate descriptor space by means of genetic algorithm selection. Table S1 reports the descriptors and parameters employed for the given GTM model.

Generative Topographic Mapping

Generative Topographic Mapping (GTM) is a dimensionality reduction method, corresponding to a probabilistic extension of Self-Organizing Maps [19] , which allows visualizing the data distribution on a 2-dimensional map. A more detailed description of GTM underlying algorithms can be found elsewhere [9] . Briefly: a squared grid of nodes is generated by inserting a flexible 2D manifold, into the initial high-dimensional space in which items occupy specific points defined by their attribute (descriptor vectors). The manifold is deformed in order to match (to approach) a maximum of these "frame" items, then it is flattened out with the above-mentioned squared grid of points, defining the latent space. The Global dataset (Figure 1, step 2) has been used to train the GTM model (i.e. to train the manifold). Genetic algorithm [20] was employed for selecting the best ISIDA descriptor space and the characteristic parameters of the GTM, as described in the next paragraph.

Once the 2D map is created, any property (here, density or class assignment) can be added as a 3rd axis forming a property landscape [9] . Here two types of landscapes are considered: (i) density landscape which assigns a color code depending on the number of compounds populating a given GTM node; (ii) class landscapes in which the color code is assigned according to the repartition of C/nC compounds.

The naming "Universal Map" refers to the GTM model showing the best overall performances for all 11 considered endpoints, see section "2.7 Ranking the performances of GTM models".

GTM's applicability domain

The "fragment control" assessment [2] is employed as method to define a model's Applicability Domain: if the test molecule has a fragment not present in the GTM's training set (i.e. the Global dataset), it is considered to be outside of applicability domain of the model. The profiling on the REACH-INV compounds has been performed only on those compounds which fulfilled the Applicability Domain requirement. Additionally, GTM models have a build-in applicability domain: any compound projected in empty regions of the chemical space are discarded. Empty regions appear as white areas on GTM landscapes (Figure 5).

Genetic algorithm optimization

A Genetic algorithm-driven optimization (Figure 1, steps 3, 4 and 5) [20] was run in order to choose the most appropriate descriptor space and GTM hyperparameters, such as the number of radial basis function centers (m), the radial basis functions width (w), the dimension of the map (k) and the regularization coefficient (l). All these parameters are encoded by a chromosome, i.e. a vector of settings needed to build a given map. The genetic algorithm therefore builds hundreds of maps based on different chromosomes. Maps' performance has been evaluated by cross-validated GTM-driven classification models for each of the considered 11 endpoints. As we are dealing with two-classes, Sensitivity (Sn), Specificity (Sp) and Balanced Accuracy (BA) parameters were computed (Table S2). The latter (BA) was chosen as scoring function for the optimization process. For each endpoint, a property-specific cross-validated BA value is returned. To obtain a more robust evaluation, this cross-validation procedure is repeated three times, and the map fitness score is based on the mean of all set-specific BAs. In the end, each map has an associated vector of 11 BA values, one per endpoint. The genetic algorithm uses the concept of the Pareto-front optimization to select the optimal set of nondominated solutions [21] and filtering redundant configurations. The procedure resulted in 28 maps with a unique set-up of descriptor spaces and GTM hyperparameters (Table S1). Some of these maps can perform well on some tasks and poorly on others. This poses two challenges: (i) select the best possible map on all tasks (paragraph 2.7) and (ii) assembling an efficient ensemble model (paragraph 2.9). To rationalize these choices, GTM models candidates went through a Rank Product [22] scoring procedure.

Ranking the performances of GTM models

The genetic algorithm run identifies a set of different manifolds, each based on a particular type of ISIDA descriptors (28 descriptor spaces were considered). Since multiple endpoints can be predicted using the same manifold, the best "all-around" map can be selected using a score measuring overall performance of the considered GTM-based classification models. In principle, a mean value for the ensemble of balanced accuracies of individual models could be used for this purpose.

A more robust score (Figure 1, step 6) can be obtained by using the "Rank Product" scoring method [22] : (i) for the given property, the considered manifolds are sorted according to their BA values; (ii) a score S, equal to the rank in the sorted list, is assigned starting from the top manifold; (iii) this process of sorting and scoring assignment is repeated for each property; (iv) the overall Rank Product is calculated as the product of each property's score (Rank Product = ∏ 𝑆 𝑖 𝑛 𝑖=1 ); (v) the map having the lowest Rank Product is selected as the best, so-called "Universal Map" (UM), reflecting its ability to have good predictive power on all considered properties. On the other hand, the wording "Optimal Map" (OM) refers to the map scored by the best BA for a given property, regardless of its performances on the others. Notice that Universal Maps result from multi-task learning procedure because all 11 endpoints were used to train the GTM manifold. In comparison, Optimal Maps result from single-task learning since only one selected endpoint was used to optimize GTM parameters.

Landscapes generation

The "best" REACH Universal Map (Figure 1, steps 7 and 8) is based on the ISIDA descriptor space IIAB(2-2) [START_REF] Ruggiu | [END_REF] , i.e. atom centered fragments with a radius of two. The following GTM hyperparameters were obtained in genetic algorithm optimization: k = 20*20; m = 7*7; w = 1.2; l = 0.02884. This map was used to visualize class landscapes and data distribution on Figures 3-5 and8-9. Class landscapes built on top 5 (out of 28 considered) manifolds were used for the REACH-INV profiling.

Consensus of GTM models

Each individual map can be used to perform predictions on all 11 properties. However, the predictive performance can be improved using a consensus model combining ensemble of individual predictors [23] . To this end, the universal maps (UMs) obtained as described in section 2.7 were combined in consensus one by one following the order defined by the Rank Product. The consensus result has been calculated by majority voting of the considered UMs. When the repartition of C/nC votes was between 40 -60 % (i.e. close to the random threshold) the result did not contribute to the performances. We observed that the performances of the consensus model are already stable after adding five maps (Figure S1).

Benchmarking of GTM with other machine learning methods

Following the strategy described in our previous works [2][3][4] , the binary consensus classification models were generated on particular training sets (one per endpoint) extracted from the Global dataset. Each consensus model is an ensemble of several individual models, based on a different descriptor spaces and/or machine learning algorithm (chosen among Random Forest, Support Vector Machine or Naïve Bayesian).

In such a way, these consensus models are optimized in terms of descriptors and methods parameters and have been trained on the same data used for GTM modelling. Therefore, 3-fold CV performances have been computed and can be directly compared with those of GTM (see section 3.6). Moreover, these models have been used to replace the missing experimental values, in order to complete the substance's environmental fate and ecotoxicological profile (see section 3.4). Table S3 report the 11 consensus models and their performances. Several other tools (such as EPI Suite or VEGA [2][3][4] ) are relevant for property prediction under REACH. However, they have not been considered for benchmarking for the following reasons:

-Not all of them can predict the compounds on all the 11 considered properties.

-Their output is not directly comparable to the GTM models, requiring additional post-processing (e.g. a continuous value needs to be categorized). The RatTox is the endpoint for which the amount of data was the highest, accounting for almost 50 % of the Global dataset. The datasets on environmental persistence (SedP, SoilP and WatP) were the smallest, with only few hundreds of compounds. There was a strong overlap of compounds between the acute aquatic toxicity datasets (AlgaeTox, DaphniaTox and FishTox): this is understandable, as for higher tonnage bands compounds shall be evaluated on all the three endpoints together to provide a complete acute aquatic toxicity evaluation. On the other hand, there was limited overlap between ready biodegradability (RB) and bioconcentration (BCF) and environmental persistence datasets. Normally, RB assays are conducted at the beginning of the registration process as, if the substance is demonstrated to be rapidly degraded in the environment, other endpoints do not need to be evaluated, and experimental testing can be therefore waived. The C/nC class repartition varies as a function of endpoint: for endocrine disruption-related properties (AR/ER binding) the number of C compounds (i.e. binders) is rather small, which is a typical situation for such biological targets. However, several fluorinated molecules are also found in a central area of the map delimited by the black dashed rectangles. This zone is populated by small molecules counting less than five atoms. • Zones 3a and 3b incorporate aliphatic and aromatic compounds mainly with the ester and ether functional groups. Zone 3a contain more aromatic molecules compared to 3b. Molecules providing both functional groups are located between these two areas. • Zones 4a and 4b agglomerate nitrogen-containing compounds. Both compounds located in 4a and 4b are characterized by the presence of nitro-containing functional groups. The structural differences between these two zones are mainly related to: (i) the presence of the diazene (N=N) moiety (absent in compounds located in 4a) and (ii) the number of hydroxyl and ketone groups (common in 4a but almost absent in 4b). Table S4 (a and b) in SI reports a list of the top-10 most similar chemicals to compounds 4a and 4b depicted in Figure 3. • The low-density regions (e.g. two molecules identified by black dots) are populated by molecules containing "rare" chemotypes which are noticeably different from other compounds from the Global set. This suggests that the applicability domain of QSARs based on public data may not include some REACH related compounds issued from the industry. This observation is consistent with our earlier studies [2,4] concerning weak performances of existing models applied to compounds of industrial context.

Results

Overview of the curated datasets

Global dataset class landscapes

Figure 5 shows the class landscapes of the Global dataset for the 11 endpoints. Blue and red areas are populated, respectively, by nC and C compounds, yellow and green colors represent mixed populated areas in which both nC and C compounds are present. For several endpoints (e.g. BCF, RB, AR/ER binding) there is a clear separation of the classes, with very few mixed areas. On the other hand, the RatTox landscape is the one that has the worst class separation, as reflected by its lower prediction performances (Table S5). Figure S3(a and b) in SI depicts an alternative version of the RatTox landscape, using the different cuf-off value of 2000 mg/kg to discriminate C/nC classes. White areas correspond to unpopulated regions which size is related to the absence of experimental data. Thus, for the series of RatTox, BCF, SoilP, WatP and SedP endpoints sorted according to the reduction of the dataset size (Table 1), the white areas on the related landscapes increase in the same order.

Ensemble of landscapes is a convenient tool of compounds profiling. As an example, Figure 5 considers two compounds: Chlordecone (CAS 143-50-0) and p-Phenylenediamine (CAS 106-50-3) depicted by star-shaped and circle-shaped black dots, respectively. In agreement with experimental data, Chlordecone is classified by the landscapes as C for 7 out of 11 endpoints (AR/ER binding, DaphniaTox, BCF, RatTox, SedP, WatP); while p-Phenylenediamine for 6 our of 11 endpoints (ER binding, AlgaeTox, DaphniaTox, FishTox, RatTox, RB). Chlordecone was used as an insecticide but was banned due to its deleterious effects on the environment, mainly related to persistence. Aniline derivates such as p-Phenylenediamine are of concern due to their acute toxicity effects on aquatic life organisms [24][25][26] .

As mentioned above, the white areas indicate the absence of experimental data for a given zone of the chemical space. Generation of new data may help to fill such empty zones and, hence, to extend the area covered by the landscape. For this purpose, in silico predictions obtained with the help of machine-learning models described in Section 2.10, have been used instead of experimental data. Notice that the predicted values outside the applicability domain of the models were discarded. This led to the decrease of the percentage of missing data over all the 11 endpoints from 89 % to 32 %. As one may see from Figure 5 (right side), the updated landscapes cover much larger area than the initial ones. The biggest improvement is observed for the smallest datasets BCF, SedP, SoilP and WatP. S5 reports cross-validation Balanced Accuracies for each Optimal Map (OM). Related GTM-based classification models demonstrates from moderate to satisfactory performances, with BAs ranging from 0.66 (RatTox) to 0.81 (BCF). The worst performing endpoints are ER binding and RatTox. The noticeable difference in BA between AR/ER binding is quite surprising, as these datasets show a high overlap of compounds. It seems that ER data are more noisy: these results are consistent with those reported by the CERAPP/CoMPARA workgroups [START_REF] Mansouri | [END_REF]13] (Table S3). For all endpoints except SedP, the sensitivity (i.e. detection of truly C compounds) is always higher than specificity. The RatTox shows the largest difference between these two metrics (Sensitivity = 0.83 and Specificity = 0.52), indicating that the RatTox map frequently misclassifies nC compared to C compounds.

Benchmarking studies.

Figure 7 and Table S6 reports cross-validated balanced accuracies for both the best UMs and for the Consensus of five selected UMs. Besides, performances were benchmarked against the 11 machine-learning models described in Section 2.10. As expected, for a given endpoint, the best Universal Maps perform worse than related Optimal Maps. On the other hand, Consensus of the top five Universal Maps prides with similar to OMs results. Compared to the models obtained with popular machine-learning methods (see Section 2.10), GTM displays slightly worse performances.

Profiling the REACH-INV dataset

The REACH-INV dataset has been profiled by the best universal map: predictions on all properties are available on Zenodo: 10.5281/zenodo.3872735. A total of 72 % of REACH-INV compounds fell inside the applicability domain of GTMbased models according to fragment control approach. Table 2 reports the number of compounds predicted as C or nC for the given property. RB is the only property for which most of the compounds were classified as C (not readily biodegradable). This was expected as ready biodegradability assays are very stringent first-tier experiments that generally underestimate the biodegradation potential. The aquatic toxicity endpoints have a similar behavior, with roughly half of the compounds predicted as C. We also found that several chemical families, such as quaternary ammonium salts, long chain alcohols and quinones were predicted toxic for all the three trophic levels (Algae, Daphnia and Fish). Only a limited amount of compounds (3-7 %) were classified of concern (C) for bioconcentration, rat toxicity and environmental persistence. To facilitate the chemical space analysis, the in-house "constrained screening" tool [27] has been used. It allows to superpose the class landscapes and, in such a way, to isolate regions of the chemical space populated by compounds possessing a given (eco)toxicological profile. Below, we provide with 3 examples of specific queries considering either all endpoints together, or only selected endpoints (e.g. aquatic toxicity or environmental fate endpoints).

Figure 9 (a,b,c) depicts the result of the superposition process, where the color code refers to the Overall Concern Score (OCS) cumulating the C labels for all considered endpoints (here, OCS varies from 0 to 11) .

In the first case (Figure 9a), all 11 landscapes shown on Figure 8 were used. No regions of the chemical space populated by the compounds labelled C with respect to endpoints were detected, the maximal OCS value was eight. Compounds located in these regions normally show acute toxicity to the aquatic environment, are not expected to rapidly degrade and in several instances exhibit acute oral toxicity. Some examples include compounds belonging the polychlorinated biphenyls (e.g. CAS 1514-82-5) which have been banned due to their deleterious effect on the environment and biota.

In the second case (Figure 9b), we focused on environmental fate landscapes (BCF and RB) aiming to extract compounds that could persist and bioconcentrate in the food chain. A consistent number of compounds belonging to the chemical family of perfluorinated compounds (e.g. CAS 118-69-4) were identified and for most of them data is scarce, especially on the bioconcentration endpoint.

In the third case (Figure 9c) the acute aquatic endpoints were considered. Chloro-and nitro-phenols (e.g., CAS 87-86-5 and 38668-48-3), some biphenyls (e.g., CAS 38668-48-3) and quaternary ammonium salts (e.g., CAS 1563-67-3) have been identified as of potential concern.

This approach is quite flexible, and it can be relevant in different regulatory contexts, for instance in order to identify Substances of Very High Concern (SVHC) under REACH. A compound is defined as SVHC if: (i) it is CMR (carcinogenic, mutagenic or toxic for reproduction); or (ii) it is PBT (persistent, bioaccumulable and toxic) or vPvB (very persistent and very bioaccumulable). In its current state, the use of reported here maps for SVHC profiling is limited, as it does not consider CMR properties nor chronic aquatic toxicity. Considered here 11 endpoints represent only a fraction of a full REACHregistration dossier, which opens a perspective to build activity landscapes for some additional endpoints. It should also be noticed that the cut-off values used to discriminate C/nC compounds should be user-controlled, in order to allow the regulator expert to adjust them depending on the context.

Conclusions

A Global dataset on 11 toxicologically-relevant endpoints resulted from the merging of multiple public data sources has been prepared. It contains: environmental fate and pathways endpoints (bioconcentration factor, ready biodegradability, environmental persistence in sediment, soil and water compartment), ecotoxicological endpoints (acute aquatic toxicity towards algae, daphnia and fish) and human health endpoints (oral acute toxicity to rats and androgen and estrogen receptor binding). A total of 17762 unique compounds listing, at least, one experimental measurement for, at least, one endpoint have been collected. Binary Concern (C) and non-Concern (nC) labels were assigned according to relevant thresholds. The Generative Topographic Mapping approach has been employed as method for data analysis and property prediction, generating single-and multitask learning models. So called REACH Universal Maps trained on the Global dataset have been generated, and their ability to classify compounds of Concern from non-Concern has been tested on all the 11 endpoints. In such a way, a given compound can be profiled on multiple properties simultaneously. The best Universal Maps display acceptable predictive performance with balanced accuracies ranging from 0.60 to 0.78, as a function of the endpoint. Assembling five best Universal Maps in a consensus improves predictive performance: balanced accuracies vary from 0.67 to 0.83. The REACH-INV dataset containing 17762 substances registered under the REACH Regulation have been profiled on the considered endpoints. Superposition of several landscapes helps to identify the zones populated by compounds of a given (eco)toxicity profile.

This work proposes a novel and unique methodology for the identification and prioritization of compounds in the context of the REACH regulation. New untested compounds can be easily profiled on several endpoints using one unique model which largely facilitates the screening process. However, as we covered only a small fraction of the properties which constitute a registration dossier. This work could be further expanded by:

-Adding more endpoints to the GTM profiler (e.g. CMR properties) to make it suitable for different regulatory frameworks, such as the identification of SVHC compounds; -Enlarging existing datasets with "waived" REACH endpoints, in particular, the outcome (C/nC) for some endpoints can be deduced from the observation of other endpoints.

-Comparing the REACH-INV chemical space versus those of other databases, such as PubChem or ChEMBL, and annotate as much as possible those databases.

-The cut-off values to discriminate C/nC compounds could be user-defined; in such a way the GTM model can be adjusted to the context (for instance, PBT or vPvB assessment).

The Global dataset and the profiled REACH-INV are available through Zenodo: 10.5281/zenodo.3872735

  [a] Laboratory of Chemoinformatics -UMR7140, University of Strasbourg, 4 Rue Blaise Pascal, 67081, Strasbourg, France [b] Toxicological and Environmental Risk Assessment unit, Solvay S.A., 85, avenue des Frères Perret, 69192, St. Fons, France

  3. REACH-INV does not contain any experimental data, but only a list of substances concerned by the Regulation. In the end, a total of 11951 compounds (out of 22966) have been retained. The Global dataset and the REACH-INV are available through Zenodo: 10.5281/zenodo.3872735.

Figure 1 .

 1 Figure 1. Modelling workflow. (1) different ISIDA descriptor spaces are generated for the Global dataset; (2) Genetic algorithm evaluates different types of descriptors and GTM parameters; (3-4) predictive performance of GTM-based models is evaluated in 3-fold crossvalidation (CV) on each of the 11 properties; (5) the Pareto front filtering is applied to exclude all dominated solutions; (6) the rank product score is applied to select the best "REACH Universal Map"; (7) the chemical space is analyzed using both density and class landscapes; (8) the REACH-INV set is profiled.

Figure 2

 2 depicts: (a) the repartition of available experimental data for 11 endpoints; and (b) the repartition of C/nC class for each particular endpoint. The Global dataset counts a total of 17762 unique compounds listing at least one experimental measurement for at least one property, for a total of 29433 data points.

Figure 2 .

 2 Figure 2. Datasets overview. (a) repartition of available experimental data for 11 endpoints; (b) Concern (C) and non-Concern (nC) classes repartition.

3. 2

 2 Density landscape: Chemical space qualitative analysis • Figure 3 depicts the density landscape of the Global dataset. The colormap emphasizes how compounds are distributed in the chemical space, identifying low-and high-densely populated regions. The chemical space is characterized according to the absolute density of datapoints in a given node, which indicates the number of compounds populating that position. The density ranges from 0 (white areas) up to 40 (yellow areas) compounds. Zone 1 is populated by aromatic and aliphatic halogenated substances, for instance belonging to the chemical family of polychlorinated biphenyl. • Zones 2a, 2b and 2c include fluorinated compounds.

Figure 3 .

 3 Figure 3. Density landscape or the Global dataset. The color scale refers to the number of compounds populating a given zone. Colored rectangles delimit map regions referred in the text. Representative structures populating selected zones of the map are shown.

3. 3

 3 Chemical space comparison: Global dataset vs. REACH-INV The inventory of substances registered under the REACHregulation (REACH-INV) has been projected on the manifold trained on the Global dataset. The likelihoods distributions of REACH-INV and the Global dataset are similar (see SI, Figure S2) which means that the manifold well describes the REACH-INV compounds. Figure 4 compares data distribution of these two databases. Blue and red colors refer to zones uniquely populated by Global dataset and REACH-INV compounds, respectively whereas intermediate colors indicate mixed regions populated by compounds from both databases. A total of 5137 out of 11951 REACH-INV compounds (43 %) are overlapping with the Global dataset. On the other hand, 12624 out of 17992 Global dataset compounds (70 %) are new to the REACH-INV. Even though the Global dataset was able to accommodate a large portion of the REACH-INV chemical space, several areas uniquely populated by REACH-INV compounds were found, indicating that the Global dataset was lacking important chemotypes: long aliphatic chain (CAS 416-630-8) and highly sulphonated compounds (in particular, belonging to the chemical class of dyes as, for instance, CAS 16470-24-9) were under-sampled in the Global dataset. The REACH-INV has also some unique chemotypes concerning perfluorinated compounds (e.g. CAS 88992-45-4).

Figure 4 .

 4 Figure 4. Global dataset and REACH-INV chemical space comparison. Blue regions are mainly populated by Global dataset compounds; red ones by the REACH-INV compounds; intermediate colors by compounds belonging to both databases. White areas display unpopulated regions.

Figure 5 .

 5 Figure 5. Class landscapes of the Global dataset for the 11 endpoints. Blue and red regions are populated by nC and C compounds, respectively. White areas correspond to unpopulated regions. Black stars and dots represent projects of two example compounds. The landscapes of the right side of the black arrows were recomputed after replacing missing experimental values with predicted ones obtained with the help of consensus classification models reported in Section 2.9. All images in a large scale are available in SI.

3. 5

 5 GTM for property prediction: single-task learning performances Figure6and Table

Figure 6 .

 6 Figure 6. Optimal Maps performances.

Figure 7 .

 7 Figure 7. Results of benchmarking studies for the 11 studied endpoints. Best UM = best-performing Universal Map; Consensus 5-UMs = the top five performing Universal Maps have been ensembled in consensus; Machine-Learning models described in Section 2.10; Optimal Map is the best performing GTM model for a given endpoint.

Figure 9 (

 9 Figure 9 (a,b,c). REACH-INV profiling maps resulted from superposition of several class landscapes shown on Figure 8. Different superposition scenarios have been considered: (a) landscapes of all 11 endpoints; (b) two environmental fate endpoints; and (c) four acute aquatic toxicity endpoints. The color code refers to the Overall Concern Score (OCS) cumulating the C labels for all considered endpoints. Some compounds identified as a concern are represented along with their CAS ID.

Table 1 .

 1 Selected endpoints and data availability.

	Endpoint	Acronym	C/nC threshold a	Unique compounds	C/nC b
	Bioconcentration factor	BCF	3 log units	1260	299 / 961

a selected threshold to discretize continuous properties into Concern (C) or non -Concern (nC) binary labels; b repartition of C and nC compounds.

Table 2 .

 2 REACH-INV GTM profiling results for the 11 endpoints predicted with the help of Universal Maps..

	Endpoint	nC	C	% (C)
	BCF	11651	300	3
	RB	3526	8425	70
	SedP	9330	2621	22
	SoilP	11581	370	3
	WatP	11077	874	7
	AlgaeTox	6972	4979	42
	DaphniaTox	6462	5489	46
	FishTox	6400	5551	47
	RatTox	11410	541	5
	AR binding	10591	1360	11
	ER binding	9984	1967	17

Figure 8. REACH-INV Class landscapes for the 11 endpoints predicted with the help of Universal Maps. All images in a large scale are available in SI.

Full Paper

C. Author et al.


Guillaume ; Bachelot
email: guillaume.bachelot@aphp.fr

Anne Bachelot

Marion Bonnier

Joe-Elie Salem

Dominique Farabos

Severine Trabado

Charlotte Dupont

Peter Kamenicky

Muriel Houang

Jean ; Fiet

Yves Le Bouc

Jacques Young

Antonin Lamazière

Combining metabolomics and machine learning

Keywords: metabolomics, 21-Hydroxylase deficiency, late onset congenital adrenal hyperplasia, hirsutism, machine learning, algorithm, Polycystic ovary syndrome, hyperandrogenism 21-Hydroxylase deficiency, late onset congenital adrenal hyperplasia, hirsutism, machine learning, algorithm, Polycystic ovary syndrome, hyperandrogenism

Participants/materials, setting, methods: Fifteen steroid species were measured in serum by liquid chromatography-mass spectrometry (LC-MS/MS). This set of 15 steroids (defined as "steroidome") used to map the steroid biosynthesis pathway was the input for our models.

Main results and the role of chance: From a single sample, modeling involving metabolic pathway mapping by profiling 15 circulating steroids allowed to perfectly identify NC21OHD from a confounding PCOS population. The constructed model using baseline LC-MS/MSacquired steroid fingerprinting successfully excluded all 17 NC21OHDs (sensitivity and specificity of 100%) from 266 PCOS from an external testing cohort of originally 549 women, without the use of ACTH testing. Blood sampling timing during the menstrual cycle phase did not impact the efficiency of our model.

Limitations, reasons for caution:

The main limitations were the use of a restricted and fully prospective cohort as well as an analytical issue, as not all laboratories are equipped with mass spectrometers able to routinely measure this panel of 15 steroids. Moreover, the robustness of our model needs to be established with a larger prospective study for definitive validation in clinical practice.

Wider implications of the findings:

This tool makes it possible to propose a new semiology for the management of hyperandrogenism. The model presents better diagnostic performances compared to the current reference strategy. The management of patients may be facilitated by limiting the use of ACTH tests. Finally, the modeling process allows a classification of steroid contributions to rationalize the biomarker approach and highlight some underlying pathophysiological mechanisms.

Study funding/competing interest(s):

This study was supported by "Agence Française de Lutte contre le dopage" and DIM Région Ile de France; This study was supported by the French institutional PHRC 2010-AOR10032 funding source and APHP). All authors declare no competing financial interests.

Introduction

21-hydroxylase deficiency (21OHD) is a rare genetic disorder with autosomal recessive inheritance caused by deleterious mutations in CYP21A2 (1). 21OHD is the most common cause (95%) of congenital adrenal hyperplasia (CAH) [START_REF] Speiser | Disease expression and molecular genotype in congenital adrenal hyperplasia due to 21-hydroxylase deficiency[END_REF]. The different clinical forms of 21OHD depend on the severity of the mutations. Besides the classical form (C21OHD), which is due to complete or severe 21OH deficiency revealed at birth, there exists the nonclassical form (NC21OHD), which is due to partial enzymatic defect, with late-onset symptoms; diagnosis is usually made in late childhood or after puberty (3). The partial preservation of cortisol and aldosterone secretion in patients with NC21OHD [START_REF] Kamenický | Cortisol and Aldosterone Responses to Hypoglycemia and Na Depletion in Women With Non-Classic 21-Hydroxylase Deficiency[END_REF] is explained by the existence of at least one of the alleles of a CYP21A2 mutation encoding a 21-hydroxylase that retains a minimum residual enzyme activity of 50% [START_REF] New | Genotype-phenotype correlation in 1,507 families with congenital adrenal hyperplasia owing to 21-hydroxylase deficiency[END_REF]. NC21OHD is most often revealed in adolescent girls and young women by a clinical pattern of hyperandrogenism accompanied by oligomenorrhea associated with increased circulating testosterone [START_REF] Admoni | Hyperandrogenism in carriers of CYP21 mutations: the role of genotype[END_REF]. Infertility and episodes of miscarriage may also occur in some women [START_REF] New | Extensive clinical experience: nonclassical 21-hydroxylase deficiency[END_REF]. The clinical presentation of patients with NC21OHD is similar to that seen in a common female condition, polycystic ovarian syndrome (PCOS), wherein signs of hyperandrogenism associated with oligomenorrhea and infertility are also very prevalent [START_REF] Matheson | Hirsutism in Women[END_REF]. As a result, many national and international consensuses recommend that the diagnosis of NC21OHD be systematically ruled out before the diagnosis of PCOS is definitively confirmed [START_REF] Teede | Recommendations from the international evidencebased guideline for the assessment and management of polycystic ovary syndrome[END_REF]. To exclude NC21OHD, many consensuses recommend that an assay of 17-hydroxyprogesterone (17OHP) be performed in the morning during the follicular phase and that the diagnosis of NC21OHD be considered if the basal serum concentration of 17OHP is increased [START_REF] Azziz | Screening for 21hydroxylase-deficient nonclassic adrenal hyperplasia among hyperandrogenic women: a prospective study[END_REF][START_REF] Speiser | A multicenter study of women with nonclassical congenital adrenal hyperplasia: relationship between genotype and phenotype[END_REF].

The thresholds proposed in the literature are variable and depend on the 17OHP assay method used-that is, either immunoassay or mass spectrometry [START_REF] Makela | Nonspecificity of a direct 17 alpha-hydroxyprogesterone radioimmunoassay kit when used with samples from neonates[END_REF]. Thus, depending on the centers and methods, the diagnosis of NC21OHD will be considered if baseline 17OHP levels are ranging between 2 and 4 ng/mL [START_REF] Azziz | Screening for 21hydroxylase-deficient nonclassic adrenal hyperplasia among hyperandrogenic women: a prospective study[END_REF][START_REF] Young | Dépistage et prise en charge du déficit en 21-hydroxylase à révélation tardive chez les femmes avec hyperandrogénie[END_REF]. In such cases, it is recommended that ACTH stimulation be performed and the diagnosis established if the 17OHP peak value is greater than 10 ng/mL [START_REF] Dewailly | Clinical and biological phenotypes in late-onset 21hydroxylase deficiency[END_REF][START_REF] Bello | Basal 17hydroxyprogesterone cannot accurately predict nonclassical congenital adrenal hyperplasia in children and adolescents[END_REF]. However, a significant percentage of women with NC21OHD may have circulating levels of 17OHP close to or even below 2 ng/mL [START_REF] Bidet | Clinical and molecular characterization of a cohort of 161 unrelated women with nonclassical congenital adrenal hyperplasia due to 21-hydroxylase deficiency and 330 family members[END_REF], 2010[START_REF] Escobar-Morreale | A Prospective Study of the Prevalence of Nonclassical Congenital Adrenal Hyperplasia among Women Presenting with Hyperandrogenic Symptoms and Signs[END_REF]. In addition, patients with increased serum 17OHP levels, but without NC21OHD, have also been reported [START_REF] Rosenfield | The Pathogenesis of Polycystic Ovary Syndrome (PCOS): The Hypothesis of PCOS as Functional Ovarian Hyperandrogenism Revisited[END_REF].

This raises the question of systematizing the ACTH test in all women with a clinical presentation compatible with PCOS syndrome, so as to exclude NC21OHD. However, medicoeconomic and practical obstacles limit the generalization of ACTH stimulation test. First, PCOS is a very prevalent disease whose frequency has been estimated at between 5 and 10% of women, whereas the prevalence of NC21OHD is much lower, estimated at nearly 1/1000 in most countries [START_REF] Kuttenn | Late-onset adrenal hyperplasia in hirsutism[END_REF]. In addition, many non-specialized centers treating women suspected of having PCOS do not have easy access to this dynamic test. On the other hand, making the diagnosis of NC21OHD is a primordial issue, as it has been demonstrated that more than 50% of patients affected by this condition carry severe loss-of-function CYP21A2 mutations in one allele [START_REF] Bidet | Clinical and molecular characterization of a cohort of 161 unrelated women with nonclassical congenital adrenal hyperplasia due to 21-hydroxylase deficiency and 330 family members[END_REF]. These patients are at high risk of giving birth to children affected by the classical form of CAH if their partner also carries a severe heterozygous mutation (probability evaluated at 1/60 people in the general population) [START_REF] Bidet | Clinical and molecular characterization of a cohort of 161 unrelated women with nonclassical congenital adrenal hyperplasia due to 21-hydroxylase deficiency and 330 family members[END_REF]. The concept of profiling or multiplexing the circulating steroidome is an emerging paradigm shift, as this method appears to be more accurate than single-steroid immunoassay quantitation to stratify or classify patients from different conditions. Specifically, in a recent preliminary work, investigators have shown that, compared to an isolated measure of serum 17OHP, simultaneous measurements of a panel of several steroids in the basal state more accurately distinguish subjects with NC21OHD from a control population (Turcu et al., 2020a).

In parallel to these advances in hormonal evaluation, machine learning (ML) approaches are emerging as interesting tools in diagnostic procedures due to their classification strength compared to classical univariate analysis traditionally used in clinical practice. As NC21OHD is an adrenal disorder and PCOS has an ovarian origin, we hypothesized that an evaluation of a steroid panel including both predominantly adrenal steroids and others with a dual adrenal and ovarian origin would help to better distinguish between these two conditions.

In this work, we combined hormonal data obtained by simultaneous measurements of a set of circulating steroids measured by liquid chromatography-mass spectrometry (LC-MS/MS) in the basal state (metabolomic approach) together withs a machine learning approach in order to build a model that could reliably distinguish women with NC21OHD from both normal women and females affected by PCOS without the use of an ACTH stimulation test.

Subjects, materials and methods

Study overview, population and data source

Study overview

As shown in Figure 1, the study was designed using a bicentric recruitment: an internal training The ethics committee ("Comité de Protection des Personnes") approved the study. All the participants signed a written informed consent.

Subjects -population

Training set (KB): Nineteen women with NC21OHD with biallelic CYP21A2 mutations were included as the training set in this study. These patients were previously evaluated in a clinical study [START_REF] Kamenický | Cortisol and Aldosterone Responses to Hypoglycemia and Na Depletion in Women With Non-Classic 21-Hydroxylase Deficiency[END_REF] performed at inclusion visit. Main clinical and genetics characteristics of these patients have been described elsewhere [START_REF] Kamenický | Cortisol and Aldosterone Responses to Hypoglycemia and Na Depletion in Women With Non-Classic 21-Hydroxylase Deficiency[END_REF]. In brief, seven of these patients carried the mild p.V281L CYP21A2 mutation in the homozygous state and 12 were compound heterozygotes with at least one minor mutation (for details see Supplementary Table SI).

Twelve patients received, at some point in their lives, hydrocortisone (10-20 mg per day) and/or dexamethasone (0.25-0.5 mg per day), prescribed transiently for important hyperandrogenism or impaired fertility, but all had received these treatments more than 1 year before inclusion.

Nineteen age-and BMI-matched heathy women, whose characteristics have been detailed

elsewhere [START_REF] Kamenický | Cortisol and Aldosterone Responses to Hypoglycemia and Na Depletion in Women With Non-Classic 21-Hydroxylase Deficiency[END_REF] [START_REF] Rosenbaum | Early central blood pressure elevation in adult patients with 21-hydroxylase deficiency[END_REF]. Exclusion criteria were: known history of cardiovascular disease and pregnancy.

Women using combined contraceptive were excluded.

All patients gave written informed consent to participate. This study was performed in accordance with the Declaration of Helsinki. Healthy volunteers of this study were recruited by the Centre d'Investigation Clinique Paris Est (CIC-1421, Pitié-Salpétrière Hospital, Paris). progesterone (P) and pregnenolone (Preg)) were measured in serum by LC-MS/MS as described elsewhere [START_REF] Fiet | A Liquid Chromatography/Tandem Mass Spectometry Profile of 16 Serum Steroids, Including 21-Deoxycortisol and 21-Deoxycorticosterone, for Management of Congenital Adrenal Hyperplasia[END_REF]. This set of 15 steroids (defined as the "steroidome") used to map the steroid biosynthesis pathway were the input for our models. All the steroid analyses were performed in the same clinical laboratory (Saint-Antoine Hospital, Paris -France). There were no missing values.

LC-MS/MS steroid profiling -data source and acquisition -variables

Data processing, machine learning algorithm implementation, performance measures and statistical analysis

After performing regular statistical analysis (mean comparison, Mann-Whitney-Wilcoxon using Python 3.7.6 environment and the SciPy Package version 1.4.1), we trained a machine learning model using the steroidome (set of 15 steroids) as input to detect (or predict)

NC21OHD from both normal and PCOS women without the use of a ACTH stimulation test.

We used orthogonal partial least squares (OPLS) as the standard metabolomic machine learning algorithm. As previously shown in Bachelot et al. and Agnani et al. (Bachelot et al., 2021;[START_REF] Agnani | A proof of concept of a machine learning algorithm to predict late-onset 21-hydroxylase deficiency in children with premature pubic hair[END_REF], we applied the OPLS algorithm for its performance on metabolomic data and the possibility to obtain easily interpretable score plot. The OPLS model was performed using SIMCA 16.0.1 software (Umetrics, Sweden).

Statistical modeling was performed using a multicentric design (two cohorts of patients from two independent clinical institutions) as described in Figure 1: the OPLS model was obtained by training the algorithm on the training set, and the performance of this model was assessed on the test set.

As a classification task, the model performances were evaluated using confusion matrix and accuracy score (the fraction of correct predictions).

Eventually, an additional evaluation was performed. We compared the performance of the model versus the standard cut-off based on the basal 17OHP concentration of 2 ng/mL and 10 ng/mL.

Results

Study cohorts -baseline characteristics

In the training set, as previously described in Kamenicky et al., controls were matched to NC21OHD for age and BMI [START_REF] Kamenický | Cortisol and Aldosterone Responses to Hypoglycemia and Na Depletion in Women With Non-Classic 21-Hydroxylase Deficiency[END_REF].

In the test set, (i.e., PSL) PCOS and NC21OHD patients exhibited similar baseline characteristics (Table I), except for oligomenorrhea, follicle number (right and left) and serum LH. As expected, these parameters were increased in PCOS compared with NC21OHD.

As a bicentric study, basal steroidomic profiles were compared in the different control (KB vs. PSL) and NC21OHD (KB vs. PSL) groups. As shown in Supplementary Table SII, levels of 17OHP, progesterone, DOC, 21-DF, corticosterone, aldosterone, and cortisone were homogeneous in patients and controls from both institutions. In the control groups, some steroid species-such as testosterone, 21 deoxycortisol, and 17OH-pregnenolone-exhibited differences between the two institutions, thereby reflecting natural inter-individual variabilities while remaining within the same orders of magnitude and within the normal value ranges. Figure 2A and figure 2B showed respectively the circulating basal concentrations of 17OHP and circulating concentrations of 17OHP after ACTH stimulation in the 2 groups of NC21OHD patients (n = 36) and PCOS patients (n = 266), which were significantly different (p < 0.0001). Some overlapping among NCAH and PCOS patients was observed for basal 17OHP concentrations. However, all PCOS (n = 266) exhibited a concentration well below the 10 ng/mL cutoff values of 17OHP after injection. Likewise, none of the NC21OHD patients (n = 36) exhibited post-stimulation 17OHP values below 10 ng/mL. Thus, there was strictly no overlap between NC21OHDs (with a minimum stimulated 17OHP concentration=16.40ng/mL) and PCOSs (with a maximum elevated 17OHP concentration=4.55ng/mL). These data are in accordance with the literature on 21-hydroxylase deficiency, which proposed that when a patient was genotyped with a biallelic mutation, the peak after ACTH was consistently above 10 ng/mL [START_REF] Bidet | Clinical and molecular characterization of a cohort of 161 unrelated women with nonclassical congenital adrenal hyperplasia due to 21-hydroxylase deficiency and 330 family members[END_REF][START_REF] Bidet | Fertility in women with nonclassical congenital adrenal hyperplasia due to 21-hydroxylase deficiency[END_REF]. These results excluded the possibility that any of the PCOS patients in the test set would be affected by 21-hydroxylase deficiency.

From model to diagnosis

The main goal of our study was to avoid the ACTH stimulation test through the use of mathematical algorithms and metabolomics. For that purpose, we built the model described in Machine learning modeling (Figure 1)

Step 1: Circulating basal steroids signatures were acquired from LC-MS/MS and preprocessing (normalization).

Step 2: OPLS-DA was performed on the training set (KB). The corresponding observation score plots (Figure 3A) confirmed the separation between NC21OHD (red dots) and control (green dots) individuals. Indeed, the horizontal direction of the score plot will capture variation between the two groups (NC21OHD and controls). The Y axis represents the decision boundary, and allowed us to define two areas: in the left quadrants, NC21OHD patients with specific steroids signature; in the right quadrants, controls with their associated steroids pattern.

Step 3: Once the OPLS-DA model had been built and trained on training set (KB), the model was then tested on an independent external test set from PSL; the corresponding individuals were plotted (black stars) regardless of their label (Figure 3B). Second, NC21OHD (red stars)

and control (green stars) were colored according to their diagnostic status (Figure 3C).

Accuracy score on independent test set was 1.00. All 17 NC21OHD and 72 controls from PSL were therefore correctly classified by the OPLS-DA model.

Extension of the algorithm to solve an additional practical clinical issue: to distinguish women

with NC21OHD from PCOS women.

As mentioned above, one important clinical issue is the differentiation between NC21OHD and polycystic ovary syndrome (PCOS) in female patients with hyperandrogenism (Table I).

Therefore, PCOS patients from PSL were tested in the model. The objective was to ensure that the algorithm was able to accurately identify and segregate the NC21OHD patients, not only from controls but also from PCOS patients. After being tested on the NC21OHD/controls model, PCOS patients were plotted on OPLS-DA (Figure 3D); all were located in the controls area (266 out 266). The Variable Importance for the Projection (VIP) plot summarized the importance (contribution) of the steroids in the model (Figure 4). The VIP plot displays the VIP values as a column plot with confidence intervals (95%) derived from jack-knifing. The most segregative basal steroid species were those with VIP values larger than 1, which included 21DF, 21DB, 11BOHA, 17OHP, A, B, and T. Interestingly, except for compound B, all these steroids precursors are upstream of the enzymatic blockade observed in subjects affected by NC21OHD. This strategy was able to set up an algorithm that is also capable of efficiently stratifying all the patients of the test set without going through an ACTH stimulation test.

Comparison of score performance of the model versus 17OHP thresholds proposed in the literature to detect NC21OHD.

To demonstrate the contribution of the model compared to the classical routine procedure, based on basal thresholds concentrations of 17OHP, we chose to compare the confusion matrix of the model against the one defined with thresholds of basal 2 ng/mL and 10 ng/mL. As shown in Table II, despite the use of LC-MS/MS (more sensitive and more specific than immunoassay), using a 17OHP cut-off value of 2 ng/mL (Table II), the test has a sensitivity of 94.1% and a specificity of 92.6% for the identification of N21HOD cases. It corresponded to one false negative and 25 false positives on the test set (PSL). Using a 17OHP cut-off value of non-stimulated 10 ng/mL (Table II), the test showed a sensitivity of 29.4% and a specificity of 100%. Based on our data, the ROC analysis (Youden's index) suggested a basal 17OHP measured by LC-MS/MS cutoff level of 1.76 ng/mL as the optimal cutoff. This cut-off provided 94.1% sensitivity, but yielded 36 false positives results, leading to a more widespread use of genetic testing.

The use of the model (Table II), exhibiting 100% sensitivity and specificity of the classical threshold, was more efficient than the use of the above-indicated thresholds.

Blood sampling timing during the menstrual cycle phase does not impact the efficiency of our model

During the luteal phase of the menstrual cycle, there is a documented increase in serum progesterone and 17OHP levels of ovarian origin [START_REF] Young | Dépistage et prise en charge du déficit en 21-hydroxylase à révélation tardive chez les femmes avec hyperandrogénie[END_REF]. These hormonal changes interfere with the diagnostic thresholds of 17OHP used. For this reason, it is recommended [START_REF] Young | Dépistage et prise en charge du déficit en 21-hydroxylase à révélation tardive chez les femmes avec hyperandrogénie[END_REF] to collect samples for the determination of these steroids during the follicular phase. However, identifying the phases of the cycle is not always easy in women with oligomenorrhea. We therefore evaluate the impact of the luteal phase of the cycle on the performance of our model.

In our set test (PSL), among 338 non-NC21OHD women, including controls and PCOS, progesterone serum levels were increased: in 31, progesterone levels were above 10 ng/mL and in 32, levels were between 5 ng/mL and 10 ng/mL, suggesting a blood sampling during the luteal phase. Twenty-five out of 31 women with increased progesterone levels also had elevated 17OHP levels (above the 2 ng/mL threshold) resulting in misclassification (false positives) using the above-indicated classic thresholds. Despite this confounding parameter, the model succeeded in correctly classifying all these women as non-affected by NC21OHD (Table II).

Discussion

To our knowledge, this proof-of-concept study is the first attempt to use a relevant steroid profile combined with a machine learning approach in order to build a model that reliably distinguishes women with NC21OHD from controls as well as PCOS women without the use of an ACTH stimulation test and regardless of the phase of the menstrual cycle.

We were confronted by a dual challenge: from a clinical point of view, NC21OHD and PCOS patients presented a similar symptomatology [START_REF] Papadakis | Polycystic Ovary Syndrome and NC-CAH: Distinct Characteristics and Common Findings. A Systematic Review[END_REF] and, from a biological/biochemical point of view, it was necessary to determine and quantify relevant biomarkers to allow a reliable identification of these patients.

A basal 17OHP level cut-off value of 2 ng/mL is recommended for screening [START_REF] Carmina | Non-classic congenital adrenal hyperplasia due to 21-hydroxylase deficiency revisited: an update with a special focus on adolescent and adult women[END_REF], although, ideally, the reference value should be established at each laboratory.

However, the election of a cut-off basal 17OHP value for the screening of 21-hydroxylase deficient in NC21OHD remains a matter of debate. Indeed, the diagnosis of NC21OHD would not have been made for 2-8% of the compound heterozygous patients carrying mild and severe CYP21A2 mutations without a stimulation test by ACTH, despite a basal 17OHP value < 2 ng/mL in large cohort [START_REF] Livadas | The spectrum of clinical, hormonal and molecular findings in 280 individuals with nonclassical congenital adrenal hyperplasia caused by mutations of the CYP21A2 gene[END_REF]. Finally, sampling timing requirements (i.e., early follicular phase of the menstrual cycle), particularly in PCOS/NC21OHD patients with irregular cycles, make the procedure even more difficult in clinical practice. In this work, the patient inclusion was designed to be as close as possible to the epidemiology of clinical hyperandrogenism. The replication cohort (i.e., the cohort used to evaluate the performance of the diagnostic strategy) was composed of 62% PCOS patients, 6% NC21OHD and 32% controls.

We have shown here that isolated basal levels of serum 17OHP cut-offs, even quantified using LC-MS/MS, do not allow sufficient sensibility and specificity to effectively distinguish NC21OHD from PCOS. Indeed, despite using the LC-MS/MS method, which is more sensitive and specific than the conventional immunological assays [START_REF] Wong | Identification of the Steroids in Neonatal Plasma that Interfere with 17 Alpha-Hydroxyprogesterone Radioimmunoassays[END_REF][START_REF] Minutti | Steroid profiling by tandem mass spectrometry improves the positive predictive value of newborn screening for congenital adrenal hyperplasia[END_REF][START_REF] Fiet | A Liquid Chromatography/Tandem Mass Spectometry Profile of 16 Serum Steroids, Including 21-Deoxycortisol and 21-Deoxycorticosterone, for Management of Congenital Adrenal Hyperplasia[END_REF], the mere application of the classic thresholds of 2 or 10 ng/mL for basal 17OHP on the patients of the replication cohort yields many misclassifications (1 false negative and 25 false positives for 2ng/mL basal 17OHP cut-off). Indeed, in the literature, studies have focused on establishing new cutoffs : among these, [START_REF] Escobar-Morreale | A Prospective Study of the Prevalence of Nonclassical Congenital Adrenal Hyperplasia among Women Presenting with Hyperandrogenic Symptoms and Signs[END_REF][START_REF] Ueland | Adrenal steroid profiling as a diagnostic tool to differentiate polycystic ovary syndrome from nonclassic congenital adrenal hyperplasia: pinpointing easy screening possibilities and normal cutoff levels using liquid chromatography tandem mass spectrometry[END_REF]. By applying these cut-offs to our patients below 2 ng/mL, false positives increased to 45 and 38 for the 1.5 and 1.7 cut-offs respectively. In accordance with these two studies, the cut-off of 1.76 based on our cohort of 355 individuals, confirmed this high rate of false positives (n=36), rendering circulating basal 17OHP as a single biomarker unsuitable for diagnosis. Furthermore, the false-negative patient remained so, independently of the chosen cut-off point. These findings strongly support the multivariate, scoring approach, which would change the semiology by using a more complete panel of the steroid biosynthetic pathway rather than a single LC-MS/MS 17OHP concentration.

The combination of a 15-steroid profile with adapted data analysis maximizes the classification performance to exclude the diagnosis of NC21OHD in the context of hyperandrogenism. Turcu et al. recently proposed a logistic regression (LR) restricted to a small set of 3 circulating steroids (17OHP, 21DF, corticosterone) measured by LCMS to discriminate between patients with NC21OHD and unaffected individuals (Turcu et al., 2020a). This preliminary monocentric study exhibited an interesting performance to distinguish healthy women from NC21OHD. However, it did not evaluate the ability of this approach to differentiate NC21OHD from a larger representative cohort of women with PCOS. In addition, the absence of an external test set did not fully allow an unbiased evaluation of the performance of the logistic regression. Considering the vast heterogeneity of PCOS clinical expression and its epidemiology, the inclusion by Turcu et al. of only 27 patients with PCOS in the group of "unaffected individuals" may have been biased and not representative of the well-recognized wide hormonal spectrum observed in women with PCOS. However, this strategy (multivariate vs 17OHP only) is mathematically relevant because it reduces the risk of spurious results (type I error/false positive and type II error/false negative) [START_REF] Wold | PLS-regression: a basic tool of chemometrics[END_REF][START_REF] Eriksson | Multi-and Megavariate Data Analysis Basic Principles and Applications[END_REF]; however, one may argue that in Turcu et al., the authors focused only on few pre-selected steroids. In our model, we used a less restricted strategy including 15 steroids to more comprehensively map the steroidal pathway, incorporating steroids that are predominantly-if not exclusively-of adrenal origin (21DB, 21DF, and 11β OH androstenedione). In our experience, these steroids are particularly adapted to compare the level hormone synthesis by each gland, distinguishing adrenal from ovarian disease [START_REF] Fiet | A Liquid Chromatography/Tandem Mass Spectometry Profile of 16 Serum Steroids, Including 21-Deoxycortisol and 21-Deoxycorticosterone, for Management of Congenital Adrenal Hyperplasia[END_REF][START_REF] Pretorius | A new dawn for androgens: Novel lessons from 11oxygenated C19 steroids[END_REF]Turcu et al., 2020b;[START_REF] Agnani | A proof of concept of a machine learning algorithm to predict late-onset 21-hydroxylase deficiency in children with premature pubic hair[END_REF]. Mathematically, a set of variables is considered to be collinear if there are one or more relations or associations between those variables. These 15 steroids are interrelated in their biosynthetic pathway, resulting in a high degree of collinearity. Consequently, and based on our previous work [START_REF] Bachelot | Proof of concept and development of a couple-based machine learning model to stratify infertile patients with idiopathic infertility[END_REF][START_REF] Agnani | A proof of concept of a machine learning algorithm to predict late-onset 21-hydroxylase deficiency in children with premature pubic hair[END_REF], we used the orthogonal partial least squares approach (OPLS): the latter is less affected by multicollinearity than logistic regression because its components are uncorrelated (i.e., orthogonal) [START_REF] Wold | PLS-regression: a basic tool of chemometrics[END_REF][START_REF] Pérez-Enciso | Prediction of clinical outcome with microarray data: a partial least squares discriminant analysis (PLS-DA) approach[END_REF][START_REF] Kjeldahl | Some common misunderstandings in chemometrics[END_REF]. On a larger set of patients (n = 366), our model allowed us to successfully exclude all (100%) patients with NC21OHD from the independent test cohort (including PCOS, controls and NC21OHD). Practically, we could reduce the dimension of our system by considering fewer variables (e.g., the four most discriminating steroids: 21DF, 21DB, 11BOHA, 17OHP) for an equally efficient NC21OHD/PCOS differential diagnosis. From a practical point of view, it must be stressed that a quantitative profiling by LC-MS/MS of the 15 steroids in a 10-minute chromatographic run is not a limiting factor and indeed facilitates its real-life clinical application [START_REF] Fiet | A Liquid Chromatography/Tandem Mass Spectometry Profile of 16 Serum Steroids, Including 21-Deoxycortisol and 21-Deoxycorticosterone, for Management of Congenital Adrenal Hyperplasia[END_REF].

Beyond the diagnostic tool, our model including 15 circulating steroids provides new insights into biochemical alterations in the context of hyperandrogenism and its debated potential ovarian/adrenal origins. The selection and ranking of variables by the algorithm are understandable given the previous results obtained on each of these molecular species, which are often considered individually. In our model, six steroids were found to be particularly significant (with VIP values >1); three metabolites are ranked with higher discriminative power than 17OHP (21DF, 21DB, and 11BOHA) and three other steroids are ranked lower (androstenedione, corticosterone, and testosterone). For years, 17OHP has been used-and is still used by many practitioners-to diagnose NC21OHD [START_REF] Carmina | Non-classic congenital adrenal hyperplasia due to 21-hydroxylase deficiency revisited: an update with a special focus on adolescent and adult women[END_REF]. However, as stated before we have shown here the interest of 21-deoxysteroid metabolites in the diagnosis of NC21OHD. Indeed, 21DF and 21DB-being strictly adrenal metabolites-may represent more specific biomarkers than 17OHP, which has a double ovarian and adrenal origin [START_REF] Turcu | Adrenal androgens and androgen precursorsdefinition, synthesis, regulation and physiologic actions[END_REF]. The exclusively adrenal origin of these steroids is likely an explanation for the significant improvement of the diagnostic performance of our model compared to the isolated use of a simple 17OHP. In women with PCOS, the main source of androstenedione is in the ovary, and this androgen is not submitted to the direct action of adrenal 11-hydroxylase; consequently, the 11β hydroxyandrostenedione serum concentration found was similar to that in normal subjects, and the ratio of androstenedione to 11-OHA was higher than in normal women [START_REF] Polson | Serum 11 β -Hydroxyandrostenedione as an Indicator of the Source of Excess Androgen Production in Women with Polycystic Ovaries*[END_REF]. Because of the ovarian and adrenal origin, basal 17OHP should be measured during the follicular phase to detect NC21OHD. Indeed, an ovarian surge of 17OHP occurs until 3 ng/mL during the luteal phase of the menstrual cycle. Moreover, slightly elevated 17OHP has been reported in women with PCOS [START_REF] Pall | The phenotype of hirsute women: a comparison of polycystic ovary syndrome and 21-hydroxylase-deficient nonclassic adrenal hyperplasia[END_REF]. The inclusion of several patients (n = 31, 9.2%) in our series who were in the luteal phase (i.e., with serum progesterone level above 10 ng/mL) suggests that our model is accurate independently of the timing of blood sampling.

Finally, in our study, using all 15 steroids, the model selected and ranked 11 betahydroxyandrostenedione as one of the most discriminating variables for isolating the NC21OHD population from the rest of the PCOS and controls without a priori knowledge.

Concentrations of 11 beta-hydroxyandrostenedione were found to be higher in NC21OHD than in PCOS and controls. The choice to measure 11OH-androstenedione to distinguish 21hydroxylase deficiency from PCOS patients appears reasonable and is indicated to add diagnostic value when corroborated with the 21 DF, 21 DB and the rest of the panel. This study has some potential limitations to consider, including the use of a restricted and fully prospective cohort. The second limitation is an analytical issue because not all laboratories, especially in France, are equipped with mass spectrometers able to routinely measure this panel of 15 steroids. Finally, the robustness of our model needs to be established with a larger prospective study for definitive validation in clinical practice. The definitive confirmation of our results will also allow us to reserve expensive and complex genetic analyses of CYP21A2 for patients with a high probability of NC21OHD. An open question that will also require further studies is the ability of our approach to distinguish NC21OHD from other causes of chronic anovulation, such as hyperprolactinemia and Cushing's syndrome, which are also differential diagnoses of PCOS [START_REF] Pugeat | Recommandations pour l'exploration des hyperandrogénies[END_REF].

The main perspective of this work would be to overcome the requirement to interrupt hormonal contraceptive treatment. In fact, the diagnosis should usually be made at least 3 months after stopping all hormonal treatment (e.g., corticosteroids, hormonal contraception); this can also be inconvenient for the patient.

Conclusion

The aim of this work was to develop a unique sampling methodology to serve as an alternative to the dynamic ACTH test to exclude the diagnosis of NC21OHD in the presence of a clinical hyperandrogenic presentation at any time of the menstrual cycle, facilitating easier management and greater patient comfort. This methodology was based on a combination of machine learning with metabolic profiling of 15 molecular species (steroids) in LC-MS/MS at the base level.

This first proof of concept is promising; however, it will require further evaluation in real clinical situations in the future. Step 1 : 15 steroid species were measured in serum by liquid chromatography-mass spectrometry (LC-MS/MS)

Figure and Table legends

Step 2 : orthogonal partial least squares (OPLS-DA) was performed on the training set (KB).

The Y axis represents the decision boundary and allowed to define 2 areas : on the left quadrants, non-classic 21OH deficiency (NC21OHD) patients with specific steroids signature, and on the right quadrants controls with their associated steroids pattern.

Step 3 : Once the OPLS-DA model built, the model was then tested on an independent external test set from PSL. After measuring their steroidome (similarly to step 1), the individuals of the test cohort were put into the model. The goal was to distinguishes women with NC21OHD from both normal and polycystic ovary syndrome (PCOS) women without the use of an ACTH stimulation test. .41 [16.06, 20.76] 19.42 [18.41, 20.43] n.s. .20 [21.07, 25.33] 22.31 [20.05, 24.57] n.s.

NC21OHD patients

Fifteen

  steroid species (21-deoxycortisol (21-DF), 11β-hydroxyandrostenedione (11β-OHA), androstenedione (D4A), testosterone (T), 17-OH-progesterone (17OHP), 17-OH-pregnenolone (17-OHPreg), dehydroepiandrosterone (DHEA), aldosterone (ALDO), cortisone (E), cortisol (F), corticosterone (B), 11-deoxycorticosterone (DOC), 21-deoxycorticosterone (21-DB),

Figure 1 ,

 1 Figure 1, which was implemented with the circulating basal steroids signature (15 steroids) and

Figure 1 :

 1 Figure 1 : Study overview, step-by-step from steroidome quantitation to explainable

Figure

  Figure 2 : Serum 17OHP levels before (A) after (B) ACTH stimulation in the two patient

Figure 3 :

 3 Figure 3 : Segregation of patients according to their basal steroid profiles.

Figure 4 :

 4 Figure 4 : Major contributing steroids to distinguish non-classic 21OH deficiency

Figure 1 :Figure 3 :Figure 4 :

 134 Figure 1 : Study overview, step-by-step from steroidome quantitation to explainable clinical decision Step 1 : 15 steroid species were measured in serum by liquid chromatography-mass spectrometry (LC-MS/MS) Step 2 : orthogonal partial least squares (OPLS-DA) was performed on the training set (KB). The Y axis represents the decision boundary and allowed to define 2 areas : on the left quadrants, non-classic 21OH deficiency (NC21OHD) patients with specific steroids signature, and on the right quadrants controls with their associated steroids pattern. Step 3 : Once the OPLS-DA model built, the model was then tested on an independent external test set from PSL. After measuring their steroidome (similarly to step 1), the individuals of the test cohort were put into the model. The goal was to distinguishes women with NC21OHD from both normal and polycystic ovary syndrome (PCOS) women without the use of a ACTH stimulation test. 335x190mm (330 x 330 DPI)

Table I : Clinical and biological characteristics of polycystic ovary syndrome (PCOS) and

 I

	stimulation, 58 patients with unknown Follicular number right, 60 patients with unknown
	Follicular number left
	§ § : 2 patients with unknown FSH, LH and Estradiol, 1 patient with unknown 17OHP one hour
	after ACTH stimulation, 6 patients with unknown Follicular number left and righ
	Data are expressed as median (interquartile range). Unpaired two-tailed Mann-Whitney U test.
	*p < 0.05; **p < 0.001; n.s. not significant
	non-classic 21OH deficiency (NC21OHD) patients from test set
	§ : 1 patient with unknown BMI, 4 patients with unknown FSH, 3 patients with unknown LH,
	2 patients with unknown Estradiol, 1 patient with unknown 17 OHP one hour after ACTH

Table II : Comparison of the discriminating power of the basal [17OHP] threshold used for screening (2ng/mL); the threshold used for diagnosis (10ng/mL) as well as the model in identifying non-classic 21OH deficiency (N21HOD) cases.

 II

	Contingency table using basal [17OHP] cut-off value of 2 ng/mL has a sensitivity of 94.1% and
	a specificity of 92.6% ; Contingency table using basal [17OHP] cut-off value of 10ng/mL
	showed a sensitivity of 29.4% and a specificity of 100%. Contingency table using the model
	developed from the 15 steroids in prime rate. The model exhibiting sensitivity and specificity
	of 100%.

Table I : Clinical and biological characteristics of polycystic ovary syndrome (PCOS) and non-classic 21OH deficiency (NC21OHD) patients from test set

 I

	§ : 1 patient with unknown BMI, 4 patients with unknown FSH, 3 patients with unknown LH, 2 patients with unknown Estradiol, 1 patient with unknown 17 OHP one hour after ACTH stimulation, 58 patients with unknown Follicular number right, 60 patients with unknown Follicular numer left NC21OHD (n=17) § § Baseline characteristics Age (years) 27.8 (23.9-31.2) 30.8 (27.8-36.9) BMI (kg/m2) 25.6 (21.3-32.4) 26 (22.4-33.7) Oligomenorrhea % 87.55 23.53 Hirsutism % 71.43 64.71 Acne % 45.49 29.41 Infertility % 15.04 17.65 Weight gain % 28.95 11.76 Baldness % 13.16 0 Follicle number right (n) 24 (19-32) 15 (9-22) Follicle number left (n) 20 (16-30) 12 (6-13) FSH (UI/L) 5.5 (4.2-6.4) 5.80 (4.0-7.2) LH (UI/L) 11.2 (7.1-16.2) 6.7 (4.4-9.4) Estradiol (pg/mL) 60.0 (45.0-100.3) 66.0 (49.0-102.3) AMH (ng/mL) 7.7 (5.2 -11.8) 5.3 (3.0 -8.1) LC-MS/MS steroids Basal 17-OH-progesterone (ng/mL) 0.60 (0.40-0.94) 6.07 (3.81-9.99) 17-OH-progesterone (ng/mL) after stimulation with ACTH 1.2 (0.9-1.6) 28.9 (25.3-45.4) Basal 21-deoxycorticosterone (ng/mL) 10.00 (10.00-10.00) 15.64 (10.00-39.45) Basal 21-deoxycortisol (ng/mL) 0.01 (0.01-0.01) 0.87 (0.37-2.11) Basal 11β-hydroxyandrostenedione (ng/mL) 1.12 (0.74-1.61) 2.64 (1.82-4.89) Basal androstenedione (ng/mL) 1.93 (1.39-2.49) 2.53 (2.04-3.50) Basal testosterone (ng/mL) 0.56 (0.41-0.74) 0.56 (0.42-0.77) Basal progesterone (ng/mL) 0.12(0.08-0.25) 0.52(0.39-3.88) Screening value : 2ng/mL [17OHP] Actual 17 NC21OHD No Yes Total No True Negative = 313 False Negative = 1 314 Predicted NC21OHD Yes False Positive = 25 True Positive = 16 41 Total 338 17 355 Diagnosis value : 10ng/mL [17OHP] Actual 17 NC21OHD No Yes Total No True Negative = 338 False Negative = 12 350 Predicted NC21OHD Yes 0 True Positive = 5 5 Total 338 17 355 Machine Learning model (15 steroids) Actual 17 NC21OHD No Yes Total No True Negative = 338 0 338 Predicted NC21OHD Yes 0 True Positive = 17 17 PCOS (n=266) § Total 338 17 355	P-value * n.s ** n.s n.s n.s n.s n.s * * n.s * n.s ** ** ** ** ** ** ** n.s **

§ § : 2 patients with unknown FSH, LH and Estradiol, 1 patient with unknown 17OHP one hour after ACTH stimulation, 6 patients with unknown Follicular number left and righ Data are expressed as median (interquartile range). Unpaired two-tailed Mann-Whitney U test. *p < 0.05; **p < 0.001; n.s. not significant

Table II : Comparison of the discriminating power of the basal [17OHP] threshold used for screening (2ng/mL); the threshold used for diagnosis (10ng/mL) as well as the model in identifying non-classic 21OH deficiency

 II 4% and a specificity of 100%. Contingency table using the model developed from the 15 steroids in prime rate. The model exhibiting sensitivity and specificity of 100%. Supplementary Table SI : Age, serum 17OHP after ACTH test and CYP21A2 genotype (allele 1 and allele 2) for each patient; Training set (KB) n = 19 and Test set (PSL) n = 17 Supplementary Table SII : basal steroidomic profiles in the different control (KB vs. PSL) and NC21OHD (KB vs. PSL) groups. (21-deoxycortisol (21-DF), 11β-hydroxyandrostenedione (11β-OHA), androstenedione (D4A), testosterone (T), 17-OH-progesterone (17OHP), 17-OH-pregnenolone (17-OHPreg), dehydroepiandrosterone (DHEA), aldosterone (ALDO), cortisone (E), cortisol (F), corticosterone (B), 11-deoxycorticosterone (DOC), 21-deoxycorticosterone (21-DB), progesterone (P) and pregnenolone (Preg)) Data are expressed as median (interquartile range). Unpaired two-tailed Mann-Whitney U test. *p < 0.05; **p < 0.001; n.s. not significant

			TRAINING SET
		Age	Serum 17OHP
	Patient #	(yrs)	1h after Synacthen®	CYP21A2 genotype
	1 2 3	36 19 29	(ng/mL) 32 25 22 Training set Allele 1 p.V281L Controls p.V281L p.V281L (n=19)	Allele 2 p.V281L p.M283V p.V281L Test set (n=72)	p-value
	4 5 17 OH progesterone 35 40 21 deoxycortisol	22 47	p.V281L p.V281L 1.00 [0.51, 1.49] 0.02 [0.01, 0.03]	p.V281L 0.93 [0.69, 1.16] p.I172N ; p.Q318X 0.01 [0.01, 0.02]	n.s. **
	6 progesterone	37	32.6	p.V281L 4.01 [1.13, 6.90]	p.V281L 3.15 [1.60, 4.70]	n.s.
	7 21 deoxycorticosterone 28	69.5	p.R356W 5.00 [4.00, 6.00]	p.G178R 5.00 [4.00, 6.00]	/
	8 11-deoxycorticosterone 38	48	p.V281L 0.04 [0.03, 0.05]	p.I172N ; p.Q318X 0.05 [0.04, 0.05]	n.s.
	9 corticosterone 10 aldosterone 11 delta4-androstenedione 40 22 30 testosterone	25 21 24.9	p.V281L 2.89 [2.01, 3.77] p.V281L 0.12 [0.07, 0.17] p.W21X 1.07 [0.85, 1.29] 0.25 [0.20, 0.29]	c.293-13C>G/A 3.94 [3.13, 4.76] 0.09 [0.07, 0.11] p.V281L 1.46 [1.27, 1.65] p.V281L 0.36 [0.31, 0.40]	n.s. n.s. * **
	12 11beta-hydroxy Delta4-34	29.7		p.V281L	p.V281L
	13 androstenedione	34	45.6	p.V281L 0.85 [0.63, 1.08]	p.P453S 1.18 [1.01, 1.36]	*
	14 dehydroepiandrosterone 34	16.4	p.V281L 4.96 [3.83, 6.10]	p.P453S 7.80 [6.41, 9.19]	*
	15 17 OH pregnenolone 22 16 32 pregnenolone 17 20 cortisol cortisone	29.8 24 32.5	p.V281L 0.97 [0.61, 1.34] p.V281L 1.87 [1.51, 2.22] p.V281L 83.69 [70.97, 96.41] 18	2.55 [1.91, 3.19] c.1-?_936+?del, p.0? 1.34 [1.13, 1.55] p.V281L 104.74 [94.34, 115.14] p.V281L	** * *
	18	21	26.5		c.293-13C>G/A	p.V281L
	19	20	74.5		p.P30L	p.R356W
			TEST SET
	Patient #	Age (yrs)	Serum 17OHP Synacthen® 1h after	CYP21A2 genotype
			(ng/mL)	Allele 1	Allele 2
	1	22	24.2		p.V281L	P30L
	2	21	41.1		p.V281L	p.V281L
	3	30	25.4		p.V281L	p.V281L
	4	28	28.3		p.V281L	p.V281L
	5	43	53.8		p.V281L	p.V281L
	6	37	41		p.V281L	p.V281L
	7	21	18.6		126C>T, 113G>A, 110T>C	126C>T, 13G>A

(N21HOD) cases. Contingency table using basal [17OHP] cut-off value of 2 ng/mL has a sensitivity of 94.1% and a specificity of 92.6% ; Contingency table using basal [17OHP] cut-off value of 10ng/mL showed a sensitivity of 29.

https://academic.oup.com/humrep Draft Manuscript Submitted to Human Reproduction for Peer Review

Draft Manuscript Submitted to Human Reproduction for Peer Review

Acknowledgements:

We thank Nicolas Helin and Bettina Ribault for their technical support.

Funding: This study was supported by Agence Française de Lutte contre le dopage and DIM Région Ile de France; This study was supported by the French institutional PHRC 2010-AOR10032 funding source and APHP).

Data availability statement : The data underlying this article will be shared on reasonable request to the corresponding author.

The ethics committee ("Comité de Protection des Personnes") approved the study.

Author's roles: G.B. performed lipidomics analysis, statistical modeling analysis, drafted the manuscript, participated in the study conception and design, interpretation of data. M.B, JE.S, D.F, S.T, C.D., P.K, M.H, J.F, Y.LB participated in acquisition of data, patients care, and in critical revision of the manuscript for intellectual content. J.Y, A.B supervised the study, participated in study conception and design, drafting the manuscript and in critical revision of the manuscript for intellectual content. A.L. supervised the study, lipidomics analysis, statistical modeling analysis, participated in study conception and design, interpretation of data and drafting the manuscript.

Conflict of interest:

All authors declare no competing financial interests.