dataesr/scientific-paragraphs-categorization huggingface dataset:
Overall, in all scRNA-seq datasets, excepting  Liu et al. (2019) , no distance captured exactly the numbers of clusters reported in the corresponding publications. However, the number of clusters inferred by OT tend to be closer to the ground-truth, while Pearson correlation tended to highly overcluster the data. Regarding the NMI and ARI scores, in two out of four datasets, OT performed better than Pearson correlation (Fig. 3B ), while for Liu scRNA and Li Tumor both OT and Pearson achieved comparable performances. In other single-cell omics data, OT still performs better than Pearson correlation in the majority of the datasets. Regarding the numbers of clusters, in Liu scATAC, both Pearson correlation and OT correctly detected the presence of three clusters. In Leukemia scATAC, composed of six cell lines, OT and Pearson correlation predicted 3 and 25 clusters, respectively.
Finally, in methylation data, both OT and Pearson underestimated the real number of clusters. As shown in Figure  3C , according to ARI and NMI values, OT showed better performances in all datasets except Leukemia data. Of note, Pearson correlation performed well according to ARI and NMI for Leukemia data, but this performance is obtained considering 25 clusters instead of the five real clusters presumably present in the data. This is a practical demonstration of the interest of inferring the optimal number of clusters based on the Silhouette score, rather than fixing it to the value reported in their original publication. The improvement provided by OT for cell-cell similarity inference (Fig. 2 ) has an impact on clustering performances. Of note, these results are in agreement with those observed when the number of clusters is fixed (Supplementary Table  S5 ).
Overall, the typical clustering workflow tends to more frequently overcluster the data (Table  1 ). Regarding ARI and NMI scores, OT shows better performances, reaching perfect score for three out of five datasets (Fig. 4A ). For scRNA-seq data, both the typical clustering workflow and OT correctly predicted the correct amount of clusters present in Liu scRNA. The two CRC datasets contain seven cell types, corresponding to our ground-truth labels. For Li Tumor, both the typical clustering workflow and OT predicted the presence of two clusters. In contrast, for Li NM, the typical clustering workflow inferred four clusters, while OT suggested the presence of six clusters. Finally, for Li cell lines, OT suggested the presence of nine clusters, while the typical clustering workflow correctly identified seven clusters.
Concerning the NMI and ARI scores, OT outperformed the typical clustering workflow in all four scRNA-seq datasets (Fig. 4B ). In other single-cell omics, OT correctly identified three clusters for Liu scATAC data, while the typical clustering workflow found seven clusters. Regarding Leukemia scATAC data, involving six cell lines, both OT and the typical clustering workflow under-clustered the data. Finally, in methylation data OT tend to be closer to the correct number of clusters. Regarding ARI and NMI (Fig. 4C ), OT always outperforms the typical clustering workflow and obtains perfect performances in Liu ATAC data. Of note, similar results are obtained when substituting the Leiden clustering with the Louvain clustering (Supplementary Text S5 and Fig. S4 ), indicating that the choice of clustering algorithm should not affect conclusions.
Open-source implementation and distribution To foster the reproducibility of all the results presented in this work, we provide a Python package and Jupyter notebook covering all the analyses performed, both available at https://github. com/ ComputationalSystemsBiology/ot-scOmics, together with all the preprocessed single-cell data used in this study. Since computing OT distances is computationally intensive, the code is designed to be run on a GPU, taking advantage of the PyTorch library. For users who do not have access to a GPU, extensive explanations are provided to run the Jupyter notebooks on the Google Collaboratory platform (http://colab. research. google. com/). Note that our code also allows for CPU computations.
Discussion In this study, we assessed the potential of OT to infer cell-cell similarities from single-cell omics data. We extensively benchmarked OT performances against state-of-the-art metrics. Interestingly, OT outperformed alternative metrics in capturing cell-cell similarities in all the 13 considered datasets. The biological relevance of this improvement was assessed by performing cell clustering. In all cases, the use of OT distance resulted in improved clustering results. Of note, different clustering algorithms have been used to test whether the observed performances were affected by the choice of the algorithm. Finally, we further challenged the metrics to detect cell-cell similarity in other single-cell omics: scATAC-seq and single-cell DNA methylation data. The improvement provided by OT is conserved also when other single-cell omics are considered, despite the high sparsity and the close-to-binary nature of scATAC-seq data (de  Souza et al.
, 2020; Xiong et al. , 2019). Single-cell datasets composed of a number of cells higher than that here employed are currently available (e. g. single-cell atlases). The scope of our analysis is to test the performances of OT and baseline measures in inferring cell-cell similarity from single-cell data. In this context, we are interested in reconstructing the distance matrices comparing all possible couples of cells present in the dataset. Such computations are quadratic in the number of cells, both in terms of memory and time complexity, thus preventing us to test performances on single-cell atlases for both baseline and OT distances. However, the computation of the complete cell-cell distance matrices is not required to cluster single-cell data, thus assuring that a clustering based on OT, or other baseline measures, can be scaled to high number of cells (e.
g. single-cell atlases). Our evaluation of metrics performances is dependent on the ground-truth labels associated with the input datasets. In cell lines, the ground-truth is well established and wide transcriptional differences exist between different cell lines. Biological samples instead consist of cell types and states having less pronounced transcriptional differences. What is a ground-truth in this case is thus less clear. We here used as ground-truth the clusters identified in the original publication of each dataset. However, such labels could be improved. In addition, given that the original publications employ Euclidean distance or Pearson correlation to define the labels, the usage of such labels as ground-truth is expected to advantage these metrics compared to alternative ones.
Acknowledgements We thank  Denis Thieffry  for the scientific feedback on the work and  Jean-Philippe Vert  for the insightful discussion during the design of this project. Funding This work was supported by the  Agence Nationale de la Recherche (ANR) - JCJC  project scMOmix and  Sanofi iTech Awards. This work was performed using HPC resources from  GENCI-IDRIS  [Grant  2021-AD011012285 ]. The work of G. Peyr e was supported by the  European Research Council  ( ERC  project  NORIA ) and the  French government  under management of  Agence Nationale de la Recherche  as part of the  'Investissements d'avenir' program , reference  ANR19-P3IA-0001  ( PRAIRIE 3IA Institute ).
Data availability Data from Liu et al. can be accessed through the supplementary data 3 and 4 of the original article or at the link https://github. com/ ComputationalSystemsBiology/momixnotebook/tree/master/data/sin gle-cell. Data from Li et al. , Corces et al. and Luo et al. can be accessed through the Gene Expression Omnibus (GEO) accession numbers: GSE81861, GSE65360 and GSE97179 respectively. The code to reproduce the analyses is available at https://github. com/ ComputationalSystemsBiology/OT-scOmics.
Updated data of 2017  [3]  indicated a decrease in VL (50,000 to 90,000 estimated cases). In France, leishmaniasis is endemic in the Mediterranean part of metropolitan France, in French Guiana, and to a lesser extent in the French West Indies. Many imported cases also are recorded. Between 1999 and 2012, the incidence of endemic leishmaniasis in the South of France was 0. 21 per 100,000 inhabitants (i. e. a mean number of 22. 6 cases per year)  [4] , and the most common notified clinical form was VL (84. 5%). Five endemic foci were described: Pyrenees-Orientales, Cevennes, Provence, Cote d'Azur, and Corsica  [5]. Infections seem to occur in two distinct environments: semi-rural hillsides with mixed forests, and the urban/ peri-urban region of Marseille  [6]. Leishmania infantum is the only species involved in the autochthonous leishmaniasis cases in France.
Data collection Cases of leishmaniasis were reported to the NRCL by more than 200 sources: university hospitals, general hospitals, health service of the French army, medical laboratories, and general practitioners. The reporting form is available at the following link: https://cnrleish. edu. umontpellier. fr/files/2018/04/Declaration_pub_2011. pdf. The collected data include: age, sex, past medical history especially causes of immunosuppression (HIV, haematological malignancy, solid organ transplant, immunosuppressive therapy, primary immune deficiency), presumed place of contamination, name of the reporting centre, clinical (e. g. number and localization of lesions) and laboratory parameters, method of diagnosis (direct examination, culture, serology, PCR), and species identification if available.
Forms are completed anonymously and are analysed in the framework of NRCL public health mission, in cooperation with Santé Publique France. Cases were defined as patients with clinical signs and symptoms compatible with at least one positive laboratory test (direct microscopic examination, culture, serology, or real-time PCR). Then, the collected cases were checked and validated by a senior parasitologist before addition to the database. For this study, data collected from 1998 to 2020 (for metropolitan France) and from 2003 to 2020 (for French Guiana and French West Indies) were used. Biological diagnosis was performed according to the protocols in use in each center: direct examination, culture, PCR (and serology for VL). Species identification was based on MLEE or molecular method or MALDI-TOF.
Acknowledgments The Working Group for the  Notification of Human Leishmanioses in France  (In alphabetical order and in addition to the main authors)  Y Balard  ( Centre National de Reference des Leishmanioses, Montpellier ),  P Delaunay  and  C Pomares  ( CHU Nice ),  S Hamane  ( CHU Saint Louis, Paris ),  H Yera  ( CHU Cochin, Paris ),  F Foulet  ( CHU de Creteil ),  S Houze ( CHU Bichat, Paris ),  X Iriart  ( CHU de Toulouse ),  A Izri  ( CHU Avicenne, Paris ),  E Lightburne  ( Services de Santedes Armees Laveran, Marseille ),  Meja Rabodonirina  ( CHU de Lyon ),  Coralie L'Ollivier  ( CHU La Timone, Marseille ),  Gloria Morizot  ( Institut Pasteur, Paris ),  Tiphaine Merguey  ( CHU de Rennes ),  Renaud Piarroux  ( CHU La Timone, Marseille ) We warmly thank all the centres that notified leishmaniasis cases to the NRCL over the years, and regret not to be able to name all our correspondents individually.
, 2011; Weng et al. , 2017) , including aggregation (Van Veen and  Kuikman, 1990; von Lützow et al. , 2007). The original framework of the Millennial model  (Abramoff et al. , 2018)  sought to update early models such as Century and Roth-C with two goals: 1) to define C pools which would be related more directly to field measurements, and 2) to reflect current understanding of soil microbial and physicochemical processes. To that end, the Millennial model included explicit representation of microbial activity, association with minerals via sorption, and aggregation of organic matter. The original version  (Abramoff et al. , 2018)  was compared with the SOM model currently used in many Earth System Models (e. g. , E3SM, CESM, ORCHIDEE), Century  (Parton et al. , 1987) , but was not tested against an independent dataset of measurements.
Parametric sensitivity, collinearity, and sensitivity to inputs For each model variant (V1, ECA, MM, LIN), we calculated a local parameter sensitivity index (S ij ), which summarizes the effect of each parameter j at timestep i for each model pool, S ij = dy i dθ j θ j y i [24] Where y is the model output and θ is the parameter value  (Soetaert and Herman, 2009). We initialized the model using the steady-state C pools at default parameter values, calculated using the stode function of R Version 4. 0. 4's rootSolve Version 1. 8. 2. 1 package  (Soetaert, 2009) , as in  Wieder et al. (2015b). Then, we ran each instance of the model for 100 years using the same repeated global average year of forcing from  Abramoff et al. (2018) , and evaluated the model output y at time step i = 100 years for each parameter value θ.
Model simulations 2. 6. 1. Data used to fit and evaluate models We evaluated the Millennial V2 and Century models using three datasets of soil fractionation measurements. For the purposes of this evaluation, we used size fractionation measurements separating MAOM from larger size fractions using a particle size upper bound that ranged from 50 to 60 μm, depending on the dataset. It is important to note that according to these size fractionation protocols, and some mixed densityand size-fractionation protocols like that of  Poeplau et al. (2017) , the MAOM fraction is defined by size alone, but the aggregate and POM fractions, though we do not isolate them in this evaluation, can be defined in different ways. Traditional fractionation protocols disperse the soil to measure POM and MAOM fractions, but measure micro-and macroaggregate fractions on soil that is minimally dispersed  (Poeplau et al.
, 2018). In most fractionation protocols, microbial biomass and LMWC may be present in both MAOM and larger size fractions, but because these pools were on average less than 1% of the total modeled SOM, we grouped them into the larger size fractions, both for simplicity and to conserve the accuracy of the MAOM fraction measurement, as this was the only fraction explicitly measured by all of the datasets. The first dataset, Viscarra Rossel (VR), is derived from 495 Australian sites described in Viscarra Rossel and Hicks (2015) and Viscarra  Rossel et al. (2019). Vegetation ranges from grass to forest with 155 grazed sites and 345 sites under minimal use or nature conservation. SOC stocks were reported in t/ha for the top 30 cm of soil and converted to g C m -2. SOC was fractionated into three fractions, originally reported as particulate organic C (POC), humic organic C (HOC), and resistant organic C (ROC), corresponding to the 50-2000 μm particle size fraction, the <50 μm particle size fraction, and charcoal measured using solid-state 13 C nuclear magnetic resonance spectroscopy  (Baldock et al.
, 2013). For the purposes of this study, we define MAOM as the operational measurement of soil particles <50 μm (HOC), while the remaining C pools are considered to be comprised of soil particles 50-2000 μm (POC). The VR dataset reported a number of soil and climate variables, including bulk density, pH, percent clay and silt, and net primary productivity (NPP) derived from the land surface model BIOS2  (Haverd et al. , 2013). Soil temperature ( • C) and soil moisture (mm 3 mm -3 ) were derived from the Global Land Data Assimilation System's (GLDAS) Noah 1 • x1 • V2. 1 land surface model output for the top 0-10 cm, 10-40 cm, and 40-100 cm of soil  (Beaudoing and Rodell, 2016). We calculated a weighted average for the three soil layers, depending on the proportion of soil in each layer which is determined by the depth of sampling at each site.
The second dataset, Georgiou (KG), is derived from 659 globallydistributed sites described in  Georgiou et al. (2021) , including forest, grass, and crops experiencing varying land management. SOC concentrations were reported in g C kg -1 soil and converted to g C m -2 using the reported depth of sampling (median: 14 cm; 1st, 3rd quartile: 10 cm, 20 cm; range: 6 cm-62 cm), and bulk density for the top 1 m of soil extracted from SoilGrids 250 m  (Hengl et al. , 2017)  aggregated to 1 km. SOC was fractionated into MAOM, defined as the soil particle size fraction <60 μm. The KG dataset reported percent clay and silt, but not pH which was derived for the top 1m of soil from SoilGrids in the same way as bulk density. NPP was derived from a 0. 5 • x0. 5 • global data product of aboveground litter production  (Li et al.
, 2019). Because NPP includes both above and belowground components, we estimated the belowground (i. e. , root) litter production from the aboveground litter production using the ratio of aboveground-to-belowground biomass carbon density (Mg/ha) from gridded global maps at 300 m spatial resolution  (Spawn et al. , 2020). This makes the simplifying assumption that the ratio of aboveground-to-belowground biomass is analogous to the ratio of aboveground litter -to-fine root litter. However, the fraction of belowground biomass in  Spawn et al. (2020)  (median: 0. 44, range: 0-1) was similar to the fraction of NPP allocated to fine roots (median: 0. 47, range: 0. 08-0. 94) in a globally-distributed dataset (N = 112)  (Xia et al. , 2019) , so we accept this first order assumption for estimating belowground litter production.
The sum of the aboveground litter production from  Li et al. (2019)  and the root production estimated from  Li et al. (2019)  and  Spawn et al. (2020)  was used as the NPP input to the Millennial V2 model for the KG dataset. Soil temperature and moisture were derived from GLDAS as described for the VR dataset. The third dataset, LUCAS, is derived from 175 European sites included in the Land Use/Cover Area Frame Survey  (Tóth et al. , 2013) , representing natural and re-vegetated forests and grassland. These sites measured soil fractions by size fractionation on topsoil with vegetation residues and litter removed, sampled from 0 to 20 cm  (Cotrufo et al. , 2019). These measurements were used to evaluate the MEMS model, another model proposing measurable soil C pools  (Robertson et al. , 2019).
The MAOM pool is defined as the <53 μm soil particle size fraction, while the remaining C pools are considered to be comprised of soil particles 53-2000 μm. The dataset also measures ancillary soil properties such as pH and percent clay and silt. Bulk density, NPP, soil temperature and soil moisture were not reported, and were derived in the same way as described for the KG dataset. All data used in this study are summarized in Table  2  and a map of site locations plotted in Fig. S1.
0. 4's FME Version 1. 3. 6. 1 package with the Levenberg-Marquardt algorithm. Parameter lower bounds were either -Inf or 0, and upper bounds were either 1 or Inf. We calculated the coefficient of determination (R 2 ) for observed and predicted pools across sites using the optimal parameter set from fitting. As a model performance metric, we calculated the root-mean-square error (RMSE) as RMSE = ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅ ∑ n i=1 ( SOC obs,i -SOC mod,i ) 2 n √ [25] Where SOC obs,i is the observed SOC at each site i, SOC mod,i is the modeled SOC at each site, and n is the number of sites. We also consider the mean absolute error (MAE) and mean bias error (MBE), MAE = ∑ n i=1 ⃒ ⃒ SOC obs,i -SOC mod,i ⃒ ⃒ n [26] MBE = ∑ n i=1 ( SOC obs,i -SOC mod,i ) n [27] We calculated the Akaike information criterion (AIC) as a performance metric, which considers not only the error but also the number of model parameters, AIC = n × ln ⎛ ⎝ ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅ ∑ n i=1 ( SOC obs,i -SOC mod,i ) 2 n √ ⎞ ⎠ + 2p [28] where p is the total number of model parameters.
For the purposes of finding the best overall parameter set and describing relationships between SOC stocks and environmental variables, we fit the model to the entirety of each dataset. This is also how we calculate the in-sample R 2 , denoted R 2 in in Tables  3  and 4. However, for the purposes of cross-validation we fit the model to a random draw of 80% of each dataset, holding back 20% for testing, repeated 10 times. Model performance indices R 2 out (out-of-sample R 2 ), RMSE, MAE, MBE, and AIC are calculated by comparing modeled values to the testing dataset and averaging across the 10 repeated fits. To classify the range of predicted and observed C stocks by biome, we    (Chen et al. , 2013). Turnover time is defined here as SOC/respiration rate, assuming that steady-state has been reached  (Sierra et al.
S4 ). The Century model relationship between SOC and NPP is linear, reflecting the first-order kinetic relationship between inputs and decomposition rates. Further, the relationship between SOC and the percentage of clay and silt increases without saturating at values of clay and silt approaching 100% (Fig. S4 ). When using all the data together (N = 1329), the Millennial V2 model predicted SOC  3 ) across sites. The Millennial V2 model also outperformed the Century model when fit to each dataset separately (Table  3 , Fig. 4 ). Of the three datasets used, both the Millennial V2 and Century models were able to explain the greatest amount of site-level variation in SOC, MAOM and other fractions for the VR dataset, and the least for the LUCAS dataset (Table  3, Table 4 ). We made very weak assumptions about the prior range of the parameters (-Inf/0, 1/Inf), yet most parameters fit to Millennial V2 were within the same order of magnitude as the default parameters and similar across datasets (Table  A1 ).
The half-saturation constant for microbial uptake had a large absolute and relative change from the default parameters after fitting. Other large absolute changes included the activation energy of uptake and decomposition, as well as the half-saturation constant for decomposition. On average, the LUCAS dataset's fitted parameters were less similar to the default parameters than those of the VR and KG datasets, suggesting that SOC stocks in the LUCAS dataset respond Collinearity is measured using the index γ with a lower bound of 1 and an upper bound of infinity, where a change in one parameter can be compensated by 1-1/γ by changing other parameters (Section 2. 4). Values >10 6 are common for models with many parameters  (Soetaert, 2016). AGG = aggregates, LMWC = low molecular weight C, MAOM = mineral-associated organic matter, MIC = microbial biomass, POM = particular organic matter.
7 , Fig. S6 ) as well as the breakdown between MAOM and non-MAOM fractions (Fig. 8 , Fig. S7 ). The Millennial V2 model predicts total SOC stocks within the range of observed values for all but the montane biome (Fig. 7 ), although there is much more variation when each dataset is considered individually (Fig. S8 ). The Century model does not estimate the C stock in each soil fraction, but Millennial V2 is generally within the range of the observed values of C stocks in the MAOM and non-MAOM fractions (Fig. 8 , Fig. S9 ). Here it is possible to identify the fraction which contributes to poor model performance in certain biomes. For example, Millennial V2 predicts MAOM C stocks well in tropical grasslands, but overestimates C stocks in the non-MAOM fraction. Conversely, Millennial V2's poor performance in the montane and desert biomes is due to the model systematically underpredicting the MAOM fraction in these biomes.