dataesr/scientific-paragraphs-categorization huggingface dataset:
Mowgli's cell embedding and clustering outperform the state-of-the-art in controlled settings derived from cell lines data. We first focused on evaluating Mowgli's embedding and clustering performance in controlled settings derived from cell lines data. To represent a panel of realistic scenarios with different distributions of cells across three groups, we applied different transformations to a simple dataset composed of three cancer cell lines profiled with scCATseq (see Figure  2A ). The scCATseq dataset provides a joint profiling of scRNA and scATAC from HCT116, HeLa-S3 and K562 cell lines  8. Unlike simulated data, this solution allows us to avoid making assumptions on the distribution of the data. Indeed, generating simulated data following a Gaussian distribution, for instance, would favor methods that approximate single-cell data with this same distribution.
For details on the generation of these datasets see Methods. We benchmarked Mowgli, Seurat v4, MOFA+ and integrative NMF based on natural metrics for embedding and clustering performance: silhouette score, Adjusted Rand Index (ARI), and purity score (see Methods). In addition, we computed UMAP visualizations for the different methods and datasets  35. As shown in Figure  2B , overall, Mowgli provides superior performance over the current state-of-the-art according to all metrics. Indeed, in all datasets except 90% sparse, Mowgli has a performance greater or equal to that of other methods. In the 90% sparse dataset, integrative NMF has a better silhouette score than Mowgli but the same ARI and purity score. These performances are confirmed by looking qualitatively at the UMAP plots in Figure  2B.
In Mixed in RNA and Mixed in both Seurat v4 confuses populations when individual omics are not sufficient to identify the three groups. Regarding dropouts, one of the most challenging features of single-cell data  46  , Mowgli shows the highest resilience with respect to the state-of-the-art. Indeed, while a sparsity of 96% is still coherent with realistic data  47  , MOFA+ and Seurat v4 confuse the three populations in the 96% sparse dataset. On the opposite, Mowgli correctly separates the three groups of cells in 96% sparse.
Mowgli's cell embeddings and clusterings are competitive with the state-of-theart in complex and heterogeneous datasets We then benchmarked Mowgli, Seurat v4, MOFA+ and integrative NMF based on their embedding and clustering performance on five paired single-cell multi-omics datasets (see Figure  3A ). Of note, these data have been already largely used to benchmark single-cell multi-omics integrative methods  15, 48, 49. The chosen datasets span different sequencing technologies, modalities, tissues and sizes: (i) Liu is a scCAT-seq cell lines dataset by Liu et al. 8  (ii) PBMC 10X is a 10X Multiome human PBMC dataset from 10X Genomics (iii) OP Multiome is a 10X Multiome human bone marrow dataset from Open Problems 50 (iv) OP CITE is a CITE-seq human bone marrow dataset from Open Problems 50 (v) BM CITE is a CITE-seq human bone marrow dataset from Stuart et al.
25. BM CITE is the larger dataset here considered, with 29,803 cells. Supplementary Table  1  lists the modalities, numbers of cells, and numbers of cell types for each dataset. For details on data preprocessing, see Methods. We benchmarked Mowgli, Seurat v4, MOFA+ and integrative NMF based on the same natural metrics used in the previous section. Since these metrics require a groundtruth annotation, we used the cell-type annotations available from the original publications of these data. In Liu the ground-truth annotations are based on the cell line of origin and thus well-defined. On the contrary, the annotations of the other datasets were computationally derived, thus biasing the evaluation toward methods closest to their annotation pipeline. For instance, the BM CITE annotation is obtained by projecting the dataset onto an atlas using Seurat v3  25.
As displayed in Figure  3B , Mowgli can handle large single-cell datasets and deliver embedding and clustering performances competitive with the state-of-the-art, especially considering the lack of absolute ground-truth annotations on most datasets here employed. In particular, according to the silhouette score, Mowgli outperforms other methods in Liu, PBMC 10X, and OP CITE. MOFA+ performs best in the other two datasets. In terms of ARI score across resolutions, Seurat v4 performs best in PBMC 10X, OP Multiome, and OP CITE. Mowgli and Seurat v4 perform comparably in the BM CITE dataset. In the Liu dataset, only MOFA+ and Mowgli reach a maximum ARI of 1. Of note, in Liu, ARIs should be compared only at low resolution, as higher resolutions lead to overclustering. In terms of purity score, Mowgli outperforms other methods in the OP CITE dataset, and it is comparable to MOFA+ in the PBMC 10X dataset.
Finally, the purity scores of all methods are comparable in the Liu dataset. The UMAP plots in Figure  3B  give a qualitative intuition of the described performance. In OP CITE, only integrative NMF and Mowgli correctly separate subpopulations of B cells (Figure  3B circled ). In BM CITE, MAIT T-cells and subpopulations of CD8+ T-cells (Figure  3B circled ) are more neatly separated in Seurat v4 than in other methods. However, as explained previously, the annotation pipeline of BM CITE might favor Seurat v4.
Mowgli improves the biological interpretability of the state-of-the-art by providing cell-type specific factors in TEA-seq data We benchmarked Mowgli with respect to MOFA+ based on its biological interpretability (see Figure  4A ). Indeed, MOFA+ is the leading single-cell multi-omics integration tool providing a user-friendly biological interpretability of its latent dimensions  14. For this benchmark we considered a TEA-seq dataset of human PBMCs, corresponding to the paired profiling of: scRNA-seq, scATAC-seq, and surface proteins  7. This dataset allows us to test the methods on more than two omic datasets, thus taking into account more complementary layers of molecular regulation. First, MOFA+ and Mowgli were independently applied for the integration of the three omics constituting the TEA-seq data.
As the dataset was not provided with an annotation of the cells, we separately clustered the embeddings obtained from Mowgli and MOFA+'s and annotated them based on gene and protein markers (see Supp Figure  1 , see Figure  4B ). We identified in this way coarse immune cell types: CD4 T-cells, CD8 T-cells, B cells, Natural Killer (NK) cells, MAIT T-cells, Monocytes and Erythroid cells. Of note, the cell type annotations obtained with the two tools agree at 97%, and match an independent RNA-based annotation obtained through Azimuth (see Supp Figure  2 ). Both methods are thus able to recover the expected cell-types through clustering of their embeddings. To then test the biological interpretability of Mowgli and MOFA+, we evaluated the specificity of the associations between their factors and the identified immune cell types.
Indeed, most recent large scale paired multi-omics data are profiled in different centers thus creating batch correction issues. In addition, Mowgli does not contain a straightforward approach to define the number of latent dimensions. This problem has been however extensively studied in NMF literature and users can rely on classical tools like the cophenetic coefficient or the elbow method. At the same time, Mowgli is fairly robust to changes in the number of latent dimensions ( ) thus suggesting that small changes in will not affect its performance. Finally, as OT is inherently expensive to compute, Mowgli requires GPU computations for the larger datasets presented in this article. However, the availability of GPUs is nowadays a standard in research centers and this will be further enhanced in the future once larger single-cell datasets will be available.
e. the rows in our dataset. In other words, denoting and two rows in our dataset, where is the mean of the elements of and is the dot product. This choice of ground cost gave best results in our previous work  31. Due to the high-dimensionality of single-cell data, we use the entropic regularization of OT, a fast and GPU-enabled approximation of classical OT computed using the Sinkhorn algorithm  74  : where the entropy is defined as. If is set to zero, (Eq3) corresponds exactly to classical OT (Eq2). Increasing values of correspond to a more diffused coupling. In previous work, we showed the entropic regularization of OT to improve similarity inference between single-cell omics profiles compared to classical notions of distance  31. As explored in  31  , entropic regularization is expected to control the systematic noise due to technical dropouts and to the stochasticity of gene expression at the singlecell level.
The parameters control the sparsity of the columns of and. In order to make these parameters more comparable across omics and datasets, we define. Default values of , , , and yielded best results across experiments. For Liu and datasets derived from Liu, we run the method with 5 factors. For other datasets, we choose 50 factors (see Supp Figure  7 ). Similarly to Rolet et al. 40  , we alternate between minimizing (Eq4) on and. One can show that these smooth minimization problems on and are equivalent to the following smooth minimization problems on new dual variables. These problems can be solved using standard optimization methods, and the method of choice is L-BFGS, a limited-memory quasi-Newton method.
• Optimizing. We solve the following smooth minimization problem: Then, we update the primal variable as follows: The column-wise softmax is defined as: • Optimizing. We solve the following smooth minimization problem: Then, we update the primal variable as follows: Here, and denote the Legendre duals of the and functions, and their smooth closed form expressions are defined in Rolet et al. 40. In the subproblems defined above, the coefficients and parametrize softmax functions, and hence control the sparsity of distributions. When tends to zero, the softmax behaves like an argmax and the distributions tend to Diracs. The code is implemented in Python and relies on PyTorch  75  for matrix operations on the GPU and on Muon  45  and Scanpy 44 to handle single-cell multimodal data.
MOFA+ We compare Mowgli to MOFA+ 14 , a variational inference method analogous to sparse PCA for multi-omic data. We use the R interface MOFA2 with default training parameters. MOFA+ provides a parameter drop_factor_threshold designed to keep only informative factors, but we found that in practice it removed important information. For example, the benchmark in Zuo and Chen  18  only kept one factor for MOFA+, which is not enough to represent cellular heterogeneity in the data. We thus choose to keep 5 factors for Liu and the datasets derived from Liu, and 30 factors for the other datasets. These parameters gave the best results overall (see Supp Figure  8 ).
Integrative NMF We implemented a baseline NMF-based integration method by concatenating the features from the different omics and solving the optimization problem with positivity constraints: We implemented this approach using the TorchNMF package. As with MOFA+, we chose the number of factors that gave the best results overall (see Supp Figure  9 ). Note that this is almost equivalent to intNMF  76  with , which minimizes instead. However, on the considered datasets, the intNMF package was too slow to be able to include it in the benchmark.
Mixed in both. We simulate a dataset where the two modalities each confuse two cell populations, but separate two others. This makes the two omics complementary. To do so we replace the RNA profiles of all HCT cells with RNA profiles of random HeLa cells. Then, we replace the ATAC profiles of all K562 cells with ATAC profiles of random HeLa cells. Sparse. We simulate high dropout noise by randomly replacing 50%, 70% ot 90% of the values with zeros. Since the data is already sparse, the final sparsity is 82%, 90% and 96%. Rare population. We simulate the presence of a rare population by keeping only 10 randomly chosen HeLa cells.
Data preprocessing All preprocessing was performed using the Scanpy 44 and Muon 45 Python packages. RNA preprocessing. Quality control filtering of cells was performed on the proportion of mitochondrial gene expression, the number of expressed genes, and the total number of counts (using Muon's filter_obs). Quality control filtering of genes was performed on the number of cells expressing the gene (using Muon's filter_var). Cells were normalized to sum to 10,000 (using Scanpy's normalize_total), then logtransformed (using Scanpy's log1p). The top 2,500 most variable genes (1,500 for the Liu dataset) were selected for downstream analysis (using Scanpy's highly_variable_genes with flavor='seurat'). ATAC preprocessing. Quality control filtering of cells was performed on the number of open peaks and the total number of counts (using Muon's filter_obs).
Quality control filtering of peaks was performed on the number of cells where the peak is open (using Muon's filter_var). In Liu, TEA, and 10X PBMC, cells were normalized to sum to 10,000 (using Scanpy's normalize_total), then log-transformed (using Scanpy's log1p). In OP Multiome, cells were normalized using TF-IDF (using Muon's tfidf) to follow the preprocessing chosen by its authors. The most variable peaks were selected for downstream analysis (using Scanpy's highly_variable_genes with flavor='seurat'). Due to differences in the data's distribution across datasets, we chose to keep 1,500 peaks in Liu, 5,000 peaks in PBMC, and 15,000 peaks in OP Multiome and TEA. ADT preprocessing. Since the number of proteins is small and the data is less noisy than RNA or ATAC, no quality control or feature selection was performed.
Data analysis Gene Set Enrichment Analysis (GSEA). The gProfiler API  77  was used through Scanpy's enrich. Custom sources GO:CC, GO:MF, GO:BP, Azimuth, and ImmuneSigDB were retrieved from the Enrichr website  78. Gene sets enriched with adjusted p-values under 0. 05 (with Bonferroni correction) were selected for further analysis. To make genes comparable, we normalized rows of the matrix to 1. The 150 top genes for every factor were then used as an unordered input to gProfiler. Motif Enrichment Analysis. Signac 79 was used to perform Motif Enrichment Analysis, using the JASPAR2022 motif database 80. To make peaks comparable, we normalized rows of the matrix to 1. The 100 top peaks for every factor were used as input to Signac's FindMotifs. The union of the top peaks across factors constitutes the background.
Visualization. To visualize the latent representation of cells in MOFA+, integrative NMF, and Mowgli's models, we computed kNN graphs (k = 20) with the euclidean distance between the cells' low-dimensional embeddings (using Scanpy's neighbors). We used these graphs to compute 2-D UMAP 35 projections (using Scanpy's umap). For Seurat v4, 2-D UMAP projections based on WNN graphs were performed using Seurat v4's function RunUMAP. Clustering. For Mowgli, integrative NMF, and MOFA+, we clustered datasets using the Leiden algorithm 36 with varying resolutions (using Scanpy's leiden). Similarly to UMAP visualization, the inputs of the Leiden algorithm were the previously computed kNN graphs. For Seurat v4, Leiden clustering was performed using Seurat v4's function FindClusters.
We used Scikit-learn's implementation adjusted_rand_score. Specificity MOFA+ was applied with 15 factors, which are enough to represent the data (see Supp Figure  8 ). Mowgli was applied with 50 factors. In both datasets, a coarse Leiden clustering was applied (using Scanpy's leiden with resolution 0. 2). In both datasets, each cluster was assigned a cell type based on the expression of the canonical gene and protein markers (see Supp Figure  1 ). To confirm this annotation, Azimuth was run on the RNA signal of the dataset (using the Azimuth web tool and the PBMC reference). The agreement of the three independent annotations is confirmed in a Sankey diagram (see Supp Figure  2 ). Dendritic cells are absent in our manual annotations because of the coarseness of the clustering. Likewise, the ADT signal (see Supp Figure  5 ) informs us that there is a CD4-CD8-T cell population missed in all three annotations.
Funding This work was supported by the  Sanofi iTech Awards. The project leading to this manuscript has received funding from the  Agence Nationale de la Recherche (ANR)   JCJC project scMOmix  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 ). The work of G. Peyré 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 ). GPU computations were performed using HPC resources from  GENCI-IDRIS  (Grant  2022-AD011012079R2 ).
Two-part joint models for a longitudinal semicontinuous biomarker and a terminal event have been recently introduced based on frequentist estimation. The biomarker distribution is decomposed into a probability of positive value and the expected value among positive values. Shared random effects can represent the association structure between the biomarker and the terminal event. The computational burden increases compared to standard joint models with a single regression model for the biomarker. In this context, the frequentist estimation implemented in the R package frailtypack can be challenging for complex models (i. e. , a large number of parameters and dimension of the random effects). As an alternative, we propose a Bayesian estimation of two-part joint models based on the Integrated Nested Laplace Approximation (INLA) algorithm to alleviate the computational burden and fit more complex models.
INTRODUCTION Estimation of joint models for longitudinal and time-to-event data was initially introduced using maximum likelihood estimation  (Wulfsohn & Tsiatis, 1997; Henderson et al. , 2000; Song et al. , 2002; Chi & Ibrahim, 2006). It was further developed within the Bayesian framework in situations where maximum likelihood estimation with asymptotic assumptions faces nonidentifiability issues. It allows flexible and more complex association structures and can handle multiple longitudinal outcomes  (Andrinopoulou & Rizopoulos, 2016). Bayesian joint models can be fitted with the R package JMbayes  (Rizopoulos, 2016) , which has been used in many biomedical researches  (Lawrence Gould et al. , 2015) , among other packages (e. g. , rstanarm;  Muth et al. , 2018). Bayesian estimation is usually based on Markov Chain Monte Carlo (MCMC) techniques  (Hanson et al.
In this article, we propose an efficient Bayesian estimation procedure for the TPJM which relies on the INLA algorithm, as implemented in the R package INLA. In practice, we used the R package INLAjoint that facilitates the fit of joint models with INLA and make it more user-friendly, as INLA was not developed specifically for this class of models. The Bayesian inference is compared to the frequentist estimation of the TPJM available in the R package frailtypack  (Król et al. , 2017). The remainder of the article is structured as follows: in Section 2, we describe the TPJM and introduce the frequentist and Bayesian estimations. In Section 3, we present a simulation study to assess the performance of these two estimation strategies in terms of bias, coverage probability, computation time, and convergence rate.
These two first simulation scenarios include 200 individuals in each dataset, corresponding to a small sample size commonly seen in randomized clinical trials while the third simulation scenario includes 500 individuals in each dataset. The last simulation scenario includes a sample size of 200 individuals and natural cubic splines for the biomarker specified with the ns() function in R, with a knot at the sample median of observed times. This generates three bases for the continuous part (one intercept and two slopes separated by the knot) and one basis for the binary part (intercept). The bases are then assumed correlated random effects in the model. For each scenario, we generate 1000 datasets; we first sample the positive longitudinal biomarker repeated measurements from a Gaussian distribution and include the zero values sampled from a binomial distribution.
The association between the probability of zero value and the positive values is given by the correlated random effects. Survival times for the terminal event are generated from an exponential baseline hazard function with a scale of 0. 2, and administrative censoring is assumed to occur at the end of the follow-up (4 years). The rate of zeros is 8% (SD = 1%) for the first three scenarios and 5% for the last scenario, which is in between what we observed in our two real datasets (12% of zeros in application 1 and 4% in application 2). A zero value observation corresponds to a patient who experienced a complete disappearance of his/her target lesions and thus is extremely informative about treatment effect. The model for data generation is given by ⎧ ⎪ ⎨ ⎪ ⎩ Logit[Prob(𝑌 ij > 0)] = 𝛼 0 + 𝑎 𝑖 + 𝛼 1 ⋅ tim𝑒 𝑗 + 𝛼 2 ⋅ tr𝑡 𝑖 + 𝛼 3 ⋅ tim𝑒 𝑗 ⋅ tr𝑡 𝑖 , E[log(𝑌 ij )|𝑌 ij > 0] = 𝛽 0 + 𝑏 0𝑖 + (𝛽 1 + 𝑏 1𝑖 ) ⋅ tim𝑒 𝑗 + 𝛽 2 ⋅ tr𝑡 𝑖 + 𝛽 3 ⋅ tim𝑒 𝑗 ⋅ tr𝑡 𝑖 + 𝜀 ij , 𝜆 𝑖 (𝑡|𝑌 ij ) = 𝜆 0 (𝑡) exp ( 𝛾 ⋅ tr𝑡 𝑖 + 𝜑 𝑎 ⋅ 𝑎 𝑖 + 𝜑 𝑏 0 ⋅ 𝑏 0𝑖 + 𝜑 𝑏 1 ⋅ 𝑏 1𝑖 ] ) , ⎡ ⎢ ⎢ ⎣ 𝑎 𝑖 𝑏 0𝑖 𝑏 1𝑖 ⎤ ⎥ ⎥ ⎦ ∼ MVN ([ 𝟎 𝟎 ] , [ 𝚺 aa 𝚺 ab 𝚺 ab 𝚺 bb ]) ≡ MVN ⎛ ⎜ ⎜ ⎜ ⎝ ⎡ ⎢ ⎢ ⎣ 0 0 0 ⎤ ⎥ ⎥ ⎦ , ⎡ ⎢ ⎢ ⎢ ⎣ 𝜎 2 𝑎 𝜎 𝑎𝑏 0 𝜎 𝑎𝑏 1 𝜎 𝑎𝑏 0 𝜎 2 𝑏 0 𝜎 𝑏 0 𝑏 1 𝜎 𝑎𝑏 1 𝜎 𝑏 0 𝑏 1 𝜎 2 𝑏 1 ⎤ ⎥ ⎥ ⎥ ⎦ ⎞ ⎟ ⎟ ⎟ ⎠.
APPLICATION We applied the Bayesian TPJM to two cancer clinical trials, the GERCOR, and the PRIME studies. A comparison with frailtypack is provided only for the GERCOR data since this approach did not converge on the PRIME study. We used the same parameterization for INLA and frailtypack, as detailed in the simulation studies but with an increased number of integration points for the Monte Carlo method in frailtypack (i. e. , 5000 points). In the context of a Bayesian approximation of the maximum likelihood estimation, we provide indications of the p-value for both frailtypack and INLA to ease the interpretation and the comparison of the results. These p-values are computed using the Z-score, assuming normal distributions. Each package has a specific criterion for model selection that takes into account the goodness of fit and the complexity of the model (i.
Description It is a randomized clinical trial investigating two treatment strategies that included a total of 220 patients with metastatic colorectal cancer. The reference strategy (arm A) corresponds to FOLFIRI (irinotecan) followed by FOLFOX6 (oxaliplatin), while arm B involves the reverse sequence. Patients were randomly assigned from December 1997 to September 1999, and the date chosen to assess overall survival was August 30, 2002. Complete data are available on 205 individuals for data analysis. Among them, 165 (80%) died during the follow-up. There are 1475 repeated measurements for the biomarker; 174 of which are zero values (12%). A summary of the dataset structure is given in Table  1. Our model uses death as the terminal event and the repeated SLD measurements (in centimeters) as the semicontinuous biomarker.