dataesr/scientific-paragraphs-categorization huggingface dataset:
4a ). Studying the geometry of these trajectories and cell distribution along them can provide exquisite insight into cell cycle speed, cell growth or even eventual cell cycle arrest. A challenging question when studying cell cycle at the single-cell level is the automatic annotation of cells with cell cycle phases. It was studied experimentally in  [24] , where the authors used genetic constructs to follow the abundance of key cell cycle proteins which they can then relate to cell cycle phases, but doing so comes with important costs and experimenter time; a natural idea would be to transfer labels from datasets annotated using this methodology to other unlabeled ones. Unfortunately, this is not as easy as it seems: differences both in preprocessing, cell types and cell cycle properties can quite drastically affect a dataset topology and geometry, making many proximity-based methods irrelevant.
A natural label transfer strategy can be pictured as follows [Fig. 4b ]. First, we carry out data integration of all datasets into a common embedding space. Then, we predict cell cycle labels of unlabeled datasets in this common space using a supervised learning approach. Finally, the learned labels can be transferred back to the original representations to be interpreted. In this experiment, we seek to automatically annotate three single-cell RNA-seq Ewing sarcoma datasets (CHLA9, CHLA10 and TC71) gathered from  [25]. To do so, we will attempt to transfer cell cycle phase, author-provided annotations contained in an osteosarcoma dataset (U2OS) gathered from  [24] , onto the three Ewing sarcoma datasets. Preprocessing differences, geometrical specificities and apparent S/G2M label mixing within the U2OS reference dataset are tough difficulties to overcome both for integration and label transfer methods [Fig.
4a ]. We first attempt to perform the integration using BKNNCorrection, setting CHLA10 as the reference dataset considering its good quality and representativity (cells are scattered uniformly around the trajectory, and the central "hole" is well resolved). Unfortunately, predicted cell cycle labels are not satisfying, almost all being predicted to belong to G1/S phases at any point of their trajectory [Sup. Fig. 4 , bottom right]. Closer inspection reveals a very small number of matches across datasets using MNN, which can be caused by an absence of orthogonality between cell cycle factors and batch effects. This is a crucial hypothesis for neighbors-based dataset integration and its lack results in a poor matching quality, making integration unreliable. Fortunately, there exist more appropriate algorithms for these situations: transportation-based matchings.
They rely on discrete optimal transport theory that have been brought into the single-cell field a few years ago in  [26] , which can be pictured as looking for the most economical way to move mass in a metric space from a point cloud onto another. This class of problems yields a natural and harmonious way to match cells across batches, by operating at the dataset level instead of operating at the item level like MNN. We can use the transmorph pre-built model TransportCorrection inspired from SCOT  [9]  and Pamona  [10] , which consists in a few preprocessing steps followed by a transport-based matching, used to project every query item onto the barycenter of its matches  [Fig 4c]. In this case we had to use the unbalanced formulation of optimal transport  [27, 28]  to account for cell cycle phase imbalance between "standard" and "fast" cell cycle datasets, which is implemented in transmorph.
Discussion Horizontal data integration and batch effect correction are key computational challenges, especially in computational biology to be able to properly analyze single-cell data from different batches or patients  [1]. We expressed the need for modular methods to tackle this problem, and demonstrated the necessity to carefully combine trustworthy cell-cell similarity algorithms with relevant embedding algorithms. We also clearly showed how deceiving data integration can be when carried out improperly, which can be extremely detrimental for subsequent analyses. This alone motivates the need for more modular tools, where every algorithmic step can be controlled if necessary. To address this need and instead of introducing another data integration technique we present transmorph, a novel modular computational framework for data integration, implemented as an open-source python library.
We provided a sound implementation for it, and demonstrated its value throughout various real-life applications both in terms of efficiency, quality and versatility. We feel that increasing modularity and user agency often leads to bloated, over-engineered and impractical pieces of software. For this reason, we provide via transmorph several pre-built integration models ready to be used in daily workflows, with high efficiency and integration quality. For more advanced and specific applications, our framework also allows building integration models from scratch by combining a variety of algorithmic modules, all of which are implemented and optimized inside our library. We eventually provide complete interfaces which allow users to implement their own computational modules if they need to. All this is endowed with a rich software ecosystem including benchmarking datasets, integration metrics, monitoring and plotting tools.
Single-cell RNA-seq datasets We use public datasets to benchmark our framework and compare its capabilities with other state-of-the-art integration pipelines. They were chosen to mimic various real-life scenarios, with total dataset sizes in the tens of thousands. All datasets contain RNA-seq data, acquired using 10X technology. • The Zhou databank was collected from  [18]  through the Curated Cancer Cell Atlas (3CA) website, and contains osteosarcoma data from 5 different patients, ranging from 866 to 14,322 cells for a total of 64,557 cells. Each cell was annotated with a cell type among chondrocyte, endothelial, fibroblast, mesenchymal stem cell (MSC), myeloid, myoblast, osteoblast, osteoclast, pericyte, T cell. • The Chen databank was collected from  [31]  using the 3CA website, and contains 61,870 nasopharyngeal cancer single-cell RNA-seq data from 14 different patients, ranging from 1,087 to 11,210 cells.
Each cell was annotated with a cell type among B cell, endothelial, epithelial, macrophage, malignant, NK cell, plasma, T cell. Raw counts have been preprocessed following standard guidelines using the scanpy python package  [15]. First, cells with low gene counts or high mitochondrial gene expression were filtered. Raw counts were then normalized to 10,000 per cell, followed by neighborhood pooling using 5 nearest neighbors. Counts were then log(1 + x) transformed, and for each dataset the top 10,000 most variable genes were kept. All these preprocessed annotated databanks can be automatically downloaded through our framework, in order to serve for benchmarking integration methods.
Cell cycle single-cell datasets and cell cycle signal Raw counts for Ewing sarcoma cell lines datasets CHLA9, CHLA10 and TC71 were obtained from  [25]. Raw counts and annotations for the osteosarcoma U2OS dataset were obtained from  [24]. They were preprocessed according to state-of-the-art guidelines. Raw counts per cell were normalized to 10,000 to account for differences in global expression, and were then log(1 + x) transformed. Top 10,000 variable genes were kept in each dataset, and U2OS/TC71 values were scaled to unit variance/zero mean, clipped to -10/10 in order to simulate differences in preprocessing between datasets coming from different sources. Data points were eventually pooled, by setting every cell counts vector to the average of its 10 nearest neighbors (neighbors were determined using euclidean distance in a 30-PC space).
We used cell cycle genes identified in  [32]  to characterize G1/S and G2/M signal. For fast cell cycle datasets TC71 and U2OS, we used only a subset of informative G1/S genes which helped to retrieve a proper loop signal (CDK1, UBE2C, TOP2A, TMPO, HJURP, RRM1, RAD51AP1, RRM2, CDC45, BLM, BRIP1, E2F8 and HIST2H2AC).
Batch k-nearest neighbors derived algorithm have been successfully applied to dataset integration in the single-cell field, for instance in the BBKNN tool  [5]. It tends to yield high quality matching results when the orthogonality of batch effect hypothesis is verified, and only needs the tuning of the k parameter (set at k = 15 in our applications). Furthermore, it also returns a much smaller number of edges compared to label matching, of the order of kna, which greatly improves performance of subsequent pipeline steps. In terms of computational efficiency, its naive algorithm needs to compute the full na × n b distance matrix which corresponds to nan b time and space complexity. This complexity being prohibitive in practice for large scale applications, we rather use nearest neighbors approximation schemes such as the nearest neighbor descent algorithm proposed in  [35] , which was found to run in quasi-linear time.
Graph embedding algorithm and clustering Joint graph embedding is an algorithm able to build a joint weighted graph whose vertices are samples from all batches, and two samples are linked together if they appear to be similar. This graph being weighted, it can then be embedded in a low dimensional space using UMAP  [16]  or MDE  [17]. The joint embedding algorithm consists of four major steps. 1. For each batch, compute its k-nn graph weighted according to UMAP membership methodology. 2. For each pair of batches, weight matching edges according to UMAP membership methodology. 3. Build a joint graph combining edges of steps 1 and 2, possibly selecting only the heaviest edges. 4. Embed the joint graph in an abstract feature space using a graph embedding optimizer such as UMAP or MDE. Let {Xi ⊂ (Xi, di), |Xi| = ni} i≤K be a set of finite datasets to integrate, each expressed in a metric space.
Let us also assume we are provided for every pair of datasets Xi and Xj a matching Mij between Xi samples and Xj samples. We start by computing, for each batch Xi its directed k-nearest neighbors graph Gi = (Xi, Ei), weighted using UMAP methodology for computing membership strength  [16]  resulting in an adjacency matrix Ki describing a k-nearest neighbors strength graph. This guarantees every sample contains at least one edge of weight 1, and helps uniformize weights regardless of batch-specific point density; note that this graph is not symmetrized yet. The next step is to convert all matching matrices to membership strength matrices so that edge weights are of the same nature as Ki k-nearest neighbors graphs. For two batches Xi and Xj associated with a matching matrix Mij ∈ R n i ×n j , we distinguish two cases.
• If (Xi, di) = (Xj, dj) which can often be achieved between datasets of the same data type using common features selection, weights can be chosen using UMAP methodology on the bipartite matching graph Mij using distance between matched points, yielding matching strength graph described by matrix Sij. • In the general case (Xi, di) ̸ = (X | , dj), we can leverage matching strength contained in Mij to produce a dissimilarity measure, for instance using the trick Dij = inv(Mij + 1) where inv(M) denotes coordinate-wise matrix inversion. This dissimilarity can then be used to compute the membership matrix as described in the previous case. Once all k-nearest neighbors strength matrices Ki and matching strength matrices Sij have been computed, they can be assembled in a joint graph G of all batches described by an adjacency matrix G whose blocks are defined as G =        K1 S12 S13.
, nK samples. Let 1 ≤ r ≤ K be the reference dataset index, and M1r ∈ R n 1 ×nr ,. , MKr ∈ R n K ×nr be K matching matrices between each batch and the reference -by convention Mrr = In r. Each batch Xs is corrected independently towards the reference Xr. Let X m s be the matched samples of Xs (Xs rows so that corresponding row in Msr contains at least one nonzero element), and X u s be the unmatched samples (the other rows). The first step is to compute the projection of each matched sample to its barycenter Y m s = F Bary Xr ,Msr (X m s ). For every matched sample, we compute each correction vector C m s = Y m s -X m s. (6) Total correction vectors Cs are then computed as follows. For every matched cell, we set the correction vector to the corresponding C m s entry. Otherwise, we set the correction vector to the one of its closest matched cell along edges of a k-nearest neighbors graph.
If there is no matched cell within the cell's connected component, we consider this may be a specific cell type with no corresponding cell in the reference dataset, and it is left uncorrected. There exist variants for this step, for instance averaging correction vectors among sets of points (e. g. neighborhood or clustering) instead of selecting just one in order to smooth the final representation. In the end, merging is simply performed as F LC Xr ,Msr (Xs) = Xs + Cs.
, sn) deviates as little as possible from the product of its marginal distributions P (s1) × P (s2),. , P (sn) where the deviation is usually defined in terms of information geometry (e. g. , as Kullback-Leibler divergence). It is shown that ICA is efficient in detecting and correcting the batch effects in omics datasets  [36]. ICA is not a data dimensionality reduction technique per se: therefore, it is usually applied on top of reduced (e. g. , by standard PCA) and whitened representation of the initial dataset. Therefore, the choice of the number of independent components is an important hyperparameter  [37]. In the simplest approach, ICA solution represents a rotation of the whitened data point cloud such that each normalized coordinate deviates as much as possible from the standard Gaussian distribution  [38].
Discrete optimal transport Discrete optimal transport (OT) problem can be naturally pictured as follows  [28]. Assuming a set of n warehouses containing goods to deliver to m factories, the optimal transport problem consists in finding the cheapest way to transport all goods to factories knowing cost of transporting goods is proportional to both mass carried and distance traveled. Originally brought into the field as a way to predict cell fate  [26] , it has more recently been shown to be an interesting asset for matching cells across datasets in integration tools like SCOT  [9]  and Pamona  [10]. Formally, let a ∈ R n and b ∈ R m be two histograms, meaning a and b coefficients are non-negative and each vector sums up to 1. In our analogy, a represents the amount of goods stored in each warehouse, and b the capacity of each factory.
We are provided a cost matrix C ∈ R n×m , where Cij is the cost of moving one unit of mass from ai to bj. The optimal transport problem from a to b given cost C can then be expressed as follows: L C (a, b) = min P∈(R + ) n×m ij CijPij s. t. P1m = a P ⊺ 1n = b (8) L C (a, b) is called the Wasserstein distance between a and b for transport cost C, and is the cheapest cost to transport all mass from a to b in this setup. The optimal transport plan P * is the valid P minimizing Eq. 8, and can be row-normalized to 1n to be used as a probabilistic matching between a and b. In practice, optimal transport can be computed between two datasets Xa ∈ R n×da and X b ∈ R m×d b with vectorized samples in row. In this case, it is common to define a = 1 n 1n and b = 1 m 1m, and to transform Xa and X b so that they are expressed in the same feature space d.
It can be achieved for instance by selecting the genes intersection, yielding Xa ∈ R n×d and Xb ∈ R m×d. The cost matrix C ∈ R n×m is then typically defined as the pairwise distance matrix between Xa and Xb , for a given distance. If Xa and X b cannot easily be embedded in a common features space, the Gromov-Wasserstein approach is in general a better alternative. We use here optimal transport as a matching algorithm, by considering the row-normalized P ij P ⊺ i 1m as the probability that cell i from dataset Xa is similar to cell j from dataset X b. There are obvious limitations to this procedure, notably the mass conservation issue: OT will always move all mass from a to b, regardless of the possible batch-specific samples. Consequently, all Xa cells will be mapped to at least one cell in X b , even though some of these Xa cells may be of a cell type missing in X b.
Even worse, if there is a class imbalance between dataset (e. g. 50% of cell type A in dataset Xa, and 25% of cell type A in dataset X b ), there will necessarily be wrong assignments using this method. Exact computation of optimal transport is furthermore computationally expensive, of the order of O((n + m) 3 ) which makes it inefficient for large scale problems (typically above 10 4 points). The next section presents an alternate formulation which provides a good approximation of the solution at a more reasonable cost, and subsampling schemes can also be used.
u and v can be approached using an iterative two-step normalization procedure. The smaller ε, the closer the objective is from unregularized formulation, at a cost of decreased convergence rate. According to  [43]  and assuming n = m for simplicity, this algorithm computes a τ -approximate solution of the original optimal transport problem in O(n 2 log(n)τ -3 ) operations. It allows tackling larger scale problems in reasonable time. In practice, the resulting transport plan is often more fuzzy and less sparse than the exact solution, which necessitates filtering small values to stay efficient. As stated previously, one of the major drawbacks of optimal transport is its constraint to always move all mass from source distribution to target distribution. As there is almost always class imbalance between single-cell datasets, this hard constraint necessarily causes matchings between cells of different cell types.
This bad property can be worked around using an alternative unbalanced optimal transport problem  [27]. The idea is to relax the hard mass conservation constraint, by rather penalizing mass discrepancy via a divergence Dφ. Given two penalty coefficients τ1 and τ2, the objective function is written as L τ C (a, b) = min P∈(R + ) n×m ij CijPij + τ1Dφ(P1m|a) + τ2Dφ(P ⊺ 1n|b) (11) This objective function can also be optimized using an adaptation of the Sinkhorn-Knopp algorithm. Removing the hard mass conservation constraint helps in practice to deal with situations of class imbalance between datasets, while staying entirely unsupervised.
Barycentric embedding and label transfer Barycentric merging is the simplest merging to set up. It works under three assumptions, (1) one batch Xr is defined as reference and all batches will be corrected towards it; (2) reference batch Xr must be expressed in a feature space; (3) for every matching Msr ∈ R ns×nr , every row must have at least one nonzero element (∥Msr1n r ∥ 0 = ns). Assumption (1) is fulfilled by user choice, and necessitates choosing a good quality batch with representative items within every sample type. Reference choice always introduces a bias in the integration, which should not be overlooked in results interpretation. Assumption (2) is easy to verify in practice, as datasets are often vectorized and represented as n × d real valued matrices. Assumption (3) necessitates to choose a semicomplete matching, which maps every sample from batch Xi to at least one sample from batch Xr.
Transportation-based matchings usually verify this assumption, while nearest neighbors-based matchings usually do not. Failing to verify assumption (3) will cause non-matched points to be projected to the 0 of Xr feature space. Let Xs be a batch to correct with respect to the reference batch Xr given a semicomplete, row-normalized matching matrix Msr. For every sample x k ∈ Xs, the k-th row α k = M sr,k provides a weighting vector which assesses the likelihood of x k corresponding to any sample of Xr. Barycentric merging F Bary will then project x k into Xr feature space Xr so that F Bary Xr ,α (x k ) = arg min x∈Xr i≤nr αi ∥Xr,i -x∥ 2 2. (12) It is easy to show x k = i≤nr αiXr,i is the solution to this problem. Therefore, F Bary can be easily generalized to project the whole Xs dataset onto Xr given Msr via F Bary Xr ,Msr (Xs) = MsrXr.
Barycentric merging has been used in several data integration pipelines such as Seurat  [8] , SCOT  [9]  and Pamona  [10] , and generally yields good results, though there are a few downside to consider. First, choosing a reference introduces a high bias in the integration, and in some applications there may be no good option for reference; for instance, every batch could miss at least one sample class. The barycenter problem also intrinsically relies on a metric. This is an issue for high dimensional problems, for instance in scRNA-seq datasets where curse of dimensionality is a real concern; in this case, barycenter has little to no interpretable sense. A common solution is to first reduce dimensionality of Xr using component analysis or non-linear methods such as UMAP  [16]  or MDE  [17].
Benchmarking methods and parameters All benchmarks have been run on a laptop equipped with 32GB of RAM, an Intel CPU i7-10750H (12 cores) processor at 5GHz and a NVIDIA GPU GeForce GTX 1650 Ti Mobile. • EmbedMNN was used on preprocessed counts with transmorph v0. 2. 0, using default parameters: "bknn" matching, 10 matching neighbors, 10 embedding neighbors, UMAP optimizer and 2 dimensions. • BKNNCorrection was used on preprocessed counts with transmorph v0. 2. 0, using default parameters: "bknn" matching, 30 matching neighbors, 10 linear correction neighbors. • TransportCorrection was used on preprocessed counts with transmorph v0. 2. 0, using solver="unbalanced", entropy epsilon=0. 02, unbalanced reg=5. • Harmony was used with default parameters directly on preprocessed counts using the rpy2 python interface.
We also tried the harmonypy python implementation, interfaced via scanpy. It successfully converged in under 10 iterations in both cases, and produced comparable results. • scVI was used on preprocessed counts following authors guidelines with n layers=2 and n latent=30, and was optimized during 124 epochs (automatically chosen by the software). • We used the scanpy implementation of BBKNN on preprocessed counts. We carried out BBKNN with default parameters on a 50-PC representation of datasets using default parameters, using neighbors within batch=3 and 10 annoy trees. • Seurat was used in RStudio after converting AnnData datasets to h5seurat using the SeuratDisk package. We carried out the integration using SelectIntegrationFeatures, FindIntegrationAnchors and Integrate-Data with default parameters.
We were not able to complete the last integration step despite our efforts due to memory usage issues. Curie, U900) for his suggestions on how to properly package and test transmorph and in alphabetical order Jonathan Bac (Institut Curie, U900), Nicolas Captier (Institut Curie, U900), Marco Ruscone (Institut Curie, U900) and Julien Vibert (former Institut Curie, U830) for the many insightful discussions and suggestions around this project.
Acknowledgements We thank  Alexander Chervov  ( Institut Curie , U900) for his suggestions on how to process TC71 and U2OS fast cell cycle datasets,  Marianyela Petrizzelli  ( Institut Curie , U900) for her help in testing transmorph on non-Linux systems,  Jane Merlevede  ( Institut Curie , U900) for testing transmorph at an early stage,  Vicent Noël  ( Institut.