CellMentor: cell-type aware dimensionality reduction for single-cell RNA-sequencing data

Domain: neuroscience
Published: 2025-12-11
Editorial outcome: major_revision

Abstract
Single-cell RNA sequencing generates high-dimensional gene expression profiles for individual cells. A key step in its analysis is dimensionality reduction, which transforms these data into lower-dimensional representations for clustering and cell type identification. Current approaches often fail to balance technical noise removal with the preservation of biologically meaningful cell-type signals. Here, we present CellMentor, a fully supervised dimensionality reduction method based on non-negative matrix factorization that integrates cell type labels directly into its optimization objective. CellMentor minimizes variation within known populations while maximizing distinctions between types, producing low-dimensional embeddings optimized for cell type identification. Across diverse simulated and experimental datasets, CellMentor shows superior cell type separation, robust batch correction, and effective detection of rare cell populations, offering a valuable tool for integrative single-cell analyses across experiments.

Introduction
Single-cell RNA sequencing (scRNA-seq) is a powerful tool for investigating cellular heterogeneity, enabling analysis of gene expression at the single-cell level. This technology has greatly enriched our understanding of biological phenomena and disease mechanisms, revealing the diversity and function of cellular systems. The established scRNA-seq analytical workflow includes quality control and filtering, data normalization, feature selection, dimensionality reduction, and unsupervised clustering1. These procedures are critical for cell type annotation, principally through the identification of marker genes concentrated in specific clusters. An alternative approach for cell-type annotation, circumventing the clustering process, employs label transfer from a reference dataset for the annotation of individual cells2. However, both strategies have intrinsic limitations. Clusters are not inherently indicative of cell types. Although unsupervised clustering can capture signals related to cell type identity, it does not necessarily ensure differentiation, particularly in the case of closely related cell types. This limitation stems from the so-called curse of dimensionality, inherent noise in scRNA-seq data due to limited mRNA capture, and gene dropouts3. Moreover, clusters frequently signify functional similarity rather than explicit cell type identity. Cell populations with similar functional properties or developmental trajectories often appear as continuous transcriptional landscapes rather than discrete clusters, challenging the fundamental assumption that clustering can effectively separate all biologically meaningful cell types. These challenges amplify when attempting to deep phenotype the data, i.e., discern subsets of similar cell types. The alternate strategy, individual cell annotation using methods like SingleR, struggles with the low sensitivity of scRNA-seq data, rendering annotation based on a restricted set of genes.

Given their high-dimensionality, noise, and sparsity, scRNA-seq datasets pose a significant challenge in analysis. Dimensionality reduction, through deconvolving gene expression into discrete latent spaces, often aids in obtaining a manageable, low-dimensional data representation. Principal component analysis (PCA) is the most common approach for dimensionality reduction in scRNA-seq analysis and serves as the default method in popular analysis packages4,5. While PCA effectively captures global variance structure, it has several limitations for scRNA-seq data: it assumes linear relationships between features, is sensitive to outliers, and importantly, cannot leverage prior biological knowledge about cell types. Furthermore, batch effects often dominate the variance in scRNA-seq data, causing PCA to emphasize technical rather than biological variation. Integration methods like Harmony6 attempt to address batch effects by aligning similar cell populations across datasets. While effective for straightforward cases, these methods often struggle with complex scenarios involving multiple batches and closely related cell types7. Deep learning approaches such as scVI8 have shown promise but require extensive computational resources and may lack interpretability. These methods, while addressing important aspects of single-cell analysis, are not specifically optimized for the cell type identification workflow that represents a core application of scRNA-seq.

Non-negative matrix factorization (NMF) has emerged as a powerful alternative. Unlike PCA, the non-negativity constraint and its ability to learn parts-based representations make it particularly well-suited for gene expression analysis9. The NMF basis vectors can reveal interpretable patterns such as co-expressed gene modules or cell-type-specific expression programs that might remain hidden when genes are analyzed individually. Several methodologies using NMF for scRNA-seq analysis have emerged, including CASSL10, which applies base NMF for semi-supervised clustering of partially labeled scRNA-seq datasets; SDCNMF11, which utilizes an NMF representation to optimize classification while imposing similarity and dissimilarity constraints; and LIGER12 uses NMF to identify shared and dataset-specific factors that facilitate integration across multiple experiments. However, these approaches use indirect supervision: LIGER uses dataset labels rather than cell type labels, CASSL applies supervision post-hoc to clustering results, and existing methods do not directly incorporate cell type discrimination into the matrix factorization objective itself. Recently, there have been advances in supervised NMF approaches. For example, the DNMF package13 implements discriminative NMF for improved classification tasks in gene expression data. Other approaches have enhanced performance by incorporating additional prior knowledge14, such as marker genes, while others have demonstrated the potential of transfer learning by using learned bases from annotated data to factorize new, unannotated datasets15.

Despite these advances, existing methods still face significant challenges in: (1) effectively balancing noise reduction with preservation of cell type-specific signals, (2) integrating data across batches without compromising biological distinctions, and (3) detecting rare cell populations.

To address these limitations, we need dimensionality reduction methods specifically designed for the cell type identification workflow—methods that create low-dimensional representations where biologically meaningful cell populations are maximally separable, facilitating accurate downstream clustering and annotation. Here, we present CellMentor, a fully-supervised NMF framework that addresses these challenges by leveraging labeled reference datasets to learn biologically meaningful latent spaces that can be transferred to new datasets. CellMentor introduces innovative approaches for rank selection along with a unique loss function that simultaneously preserves cell type identity and reduces technical variation. Our method leverages labeled reference datasets to learn biologically meaningful latent spaces that can be effectively transferred to new datasets. We demonstrate that CellMentor effectively captures fundamental biological differences between cell types while minimizing technical artifacts, enabling more accurate downstream analyses, including cell type identification.

Results
Overview of CellMentor
CellMentor is a fully supervised learning framework that directly incorporates cell type labels into the NMF optimization objective, distinct from existing semi-supervised approaches that apply supervision indirectly or post-hoc to unsupervised factorization results. The framework operates in two phases (Fig. 1). In the decomposition phase, CellMentor learns cell-type-specific patterns from a labeled reference dataset using supervised NMF. This process yields two matrices: W (genes × factors) representing gene combinations or meta-genes, and H (factors × cells) containing the cell-specific usage of these factors.

For the decomposition phase, we developed a supervised framework that incorporates discriminative constraints to preserve cell type identity. The supervised framework incorporates an optimization function with key parameters that balance reconstruction accuracy, sparsity, orthogonality, within-class compactness, and between-class separation.

In the projection phase, CellMentor applies these learned patterns to analyze new datasets. Specifically, we use the W matrix learned during decomposition to project query data into the learned latent space. This transfer learning approach enables robust cell type identification across experimental batches without requiring retraining, making it particularly valuable for integrating new datasets with established reference atlases. The resulting reduced-dimensional representation effectively manages technical variation, facilitating downstream analyses such as cell type identification and characterization of rare populations.

Development and optimization of the decomposition phase
A critical challenge in dimensionality reduction is determining the optimal number of dimensions (factors, K) that capture the true biological signal while excluding technical noise. In the context of NMF, this rank determines the number of meta-genes (columns) in the W matrix and their corresponding usage patterns (rows) in the H matrix. CellMentor addresses this through an approach combining the biwhitening method with eigenvector analysis16. The biwhitening method transforms the data matrix such that both its rows and columns have unit variance, enabling more accurate estimation of the matrix rank.

Our approach extends the biwhitening framework by incorporating two key components. First, we evaluate different scaling techniques (Poisson and quadratic variance) before applying biwhitening, selecting the approach that minimizes the distance between the empirical eigenvalue distribution and the theoretical Marchenko–Pastur (MP) distribution. The quadratic variance scaling technique consistently showed superior performance in our simulations. Second, we refine the rank estimation through eigenvector localization analysis. Using simulated data where ground truth is known, we showed that delocalized eigenvectors often correspond to technical variation rather than biological signal (Supplementary Fig. 1). By setting a significance level of 0.01 for eigenvector localization, we effectively filtered out dimensions dominated by technical noise while preserving biologically meaningful variation.

The core innovation of CellMentor lies in its supervised learning framework, which guides dimensionality reduction to explicitly preserve cell type identity. Our approach incorporates discriminative constraints into the NMF optimization, designed to maximize separation between different cell types while maintaining coherence within cell populations. This is achieved through an optimization function that balances four key components: (1) a reconstruction term that ensures the factorization accurately represents the original data, (2) a within-class compactness term that minimizes variation among cells of the same type (3) a between-class separation term that maximizes distances between different cell types, with stronger weights for similar cell types that are more likely to be confused, and (4) a sparsity constraint that promotes sparse factor representations. Each component is controlled by tunable parameters that adapt to local and global data characteristics.

Evaluation of CellMentor in simulated data
We first evaluated CellMentor using a systematic simulation framework designed to test its performance under controlled conditions with increasing levels of complexity. Using Splatter17, we generated two sets of simulated scRNA-seq datasets with known ground truth cell types and batch effects. Our objective is to evaluate the effectiveness of the different methods to create a reduced dimension for clustering and, in turn, for cell type identification. Thus, we assess performance using clustering metrics, including Adjusted Rand Index (ARI) and Normalized Mutual Information (NMI), as accurate clustering in reduced-dimensional space is the key prerequisite for reliable cell type annotation in single-cell workflows. All performance evaluations are based on clustering performed in the full reduced-dimensional space learned by each method, with UMAP visualizations provided for visualization purposes only.

In our first simulation experiment, we generated 10 distinct scRNA-seq datasets of increasing difficulty, each containing two batches and six cell types with varying degrees of technical noise. For each simulation, we designated one batch as reference data for the decomposition phase and one as query data for the projection phase, mimicking real-world scenarios where reference and query data come from different experimental conditions. We compared CellMentor against a comprehensive set of established dimensionality reduction methods spanning multiple supervision paradigms: unsupervised approaches including PCA (the standard approach used in Seurat4) and several NMF and PCA variants (ICA18, NMF19, Probabilistic Count Matrix Factorization (pCMF)20, GLM-PCA21); semi-supervised approaches specifically designed for single-cell data (CASSL10, rssNMF14, LIGER12, scCoGAPS22, scANVI8); and, for completeness, a supervised discriminant analysis method originally developed for multiview bulk data, SIDA23, adapted here to the single-cell setting. Each simulation was repeated five times with different random seeds to assess statistical robustness and method stability.

The reduced-dimensional representations obtained from different methods showed varying degrees of success in preserving cell type structure (Fig. 2A and Supplementary Fig. 2). CellMentor achieved superior performance across all simulation scenarios, with a mean ARI of 0.96 (range: 0.73–1.0), substantially outperforming other supervised methods such as CASSL (mean ARI of 0.89) and unsupervised approaches like PCA/Seurat (mean ARI of 0.719) and scANVI (mean ARI of 0.57, range: 0.48–0.803). Importantly, CellMentor demonstrated exceptional robustness to increasing simulation difficulty. While competing methods showed dramatic performance drops in simulations 7–8 and most achieved ARI scores below 0.5 in the most challenging scenarios (simulation 10), CellMentor maintained consistently high performance across all difficulty levels. In addition, some methods exhibited extreme variability between random seeds, achieving near-perfect separation in certain runs but complete failure in others (i.e., SIDA). This superior stability, combined with minimal variability across repeated runs, demonstrates CellMentor’s robustness to technical noise and biological complexity.

In our second simulation experiment, we increased the complexity by generating datasets with four batches instead of two, using two batches for reference and two for query (Fig. 2B and Supplementary Fig. 3). This scenario represents the common challenge of integrating data across multiple experimental batches while preserving cell type identity, and directly tests whether supervised dimensionality reduction might erroneously incorporate batch differences when maximizing cell type separation. We included Harmony6, a popular batch correction method, and the Siamese-based domain adaptation variant (SiameseDA)24 in our comparisons for this more complex scenario. The multi-batch scenario proved substantially more challenging for all methods, with most showing increased variability and declining performance compared to the single-batch experiments. Among supervised methods, CASSL again maintained reasonable performance in simpler scenarios but showed substantial degradation in complex simulations, while rssNMF demonstrated high variability across repetitions. Standard unsupervised methods like PCA/Seurat performed poorly throughout, with a mean ARI of only 0.67 (range: 0.57–0.72) across all simulations. Even methods specifically designed for batch integration, such as Harmony, showed inconsistent performance with a mean ARI of 0.78 and declining to 0.16 in the most challenging scenarios. CellMentor demonstrated remarkably stable and superior performance in this multi-batch scenario, maintaining an ARI of 1 across all simulations. While our simulations suggest that the projection phase may help mitigate some batch effects present in the reference data (Supplementary Fig. 4), this observation requires further validation across diverse experimental scenarios, as the effectiveness likely depends on the specific characteristics of both reference and query datasets.

To further evaluate CellMentor’s capabilities and limitations, we conducted additional simulation experiments exploring challenging scenarios in single-cell analysis (Fig. 3A). First, we tested its behavior when encountering novel cell types by excluding one cell type from the referencedataset while maintaining it in the query data (Fig. 3B). While CellMentor did not create a completely distinct cluster for this novel cell type (which would be expected for a supervised approach with no prior knowledge of this class), it still showed some ability to separate these cells from known cell types in the query data, demonstrating partial robustness to incomplete reference data.

In a related experiment, we tested how CellMentor would handle cases where reference data annotation is not sufficiently granular—a common challenge in single-cell analysis (Fig. 3C). This addresses the fundamental concern that minimizing within-class scatter might obscure biologically relevant heterogeneity. We labeled cells of one type as another type in the referencedataset, effectively merging two distinct populations under a single label. Despite training with this corrupted referencedata, CellMentor still managed to partially separate these cell populations in the query data, demonstrating that the method can detect underlying biological variation even when supervised labels are incomplete or insufficiently granular.

To assess CellMentor’s robustness to reference annotation quality, we systematically introduced mislabeling using two complementary approaches that simulate different types of real-world annotation errors. First, we tested random mislabeling by randomly reassigning 5%, 10%, 20%, 30%, 40%, and 50% of reference cells to incorrect labels (Fig. 3D, E and Supplementary Fig. 5). This experiment simulates scenarios where reference datasets contain random annotation errors, a common challenge in single-cell analysis. As expected for a supervised method, CellMentor’s performance showed gradual degradation with increasing mislabeling rates, with noticeable decreases beginning at 30% mislabeling. However, the method demonstrated remarkable robustness, maintaining an ARI score of 0.9 even with 50% of reference cells mislabeled.

To further test robustness under more biologically realistic conditions, we implemented a correlated mislabeling framework in which boundary cells close to other populations were reassigned and their expression profiles blended with the centroid of the mislabeled type (Fig. 3F). This creates transcriptional ambiguity that mimics real-world annotation errors where cells near decision boundaries are most likely to be misclassified, rather than random label swapping. Under this more challenging scenario, CellMentor still demonstrated substantial robustness: performance decreased gradually with increasing mislabeling but remained high across random seeds, with ARI values around 0.8 even at 50% mislabeling and 0.75 at 90% mislabeling. Together, these experiments demonstrate that the method can extract meaningful biological signals even from noisy annotations, whether arising from random errors or systematic misclassification of transcriptionally ambiguous cells.

We next evaluated CellMentor under conditions where the very features used for discrimination are systematically corrupted. To simulate targeted perturbations, we identified the top discriminatory genes for each cell type and progressively applied batch-specific shifts to these features, varying both the fraction of affected cell types and the fraction of their discriminatory genes. Results across the full parameter grid (Fig. 3H) show that CellMentor maintains consistently high accuracy, with mean ARI values above 0.8 even when all discriminatory features of all cell types were perturbed. In contrast, PCA and Harmony degraded sharply under the same conditions, with ARI values falling below 0.2 and 0.5, respectively. UMAP embeddings of the most challenging setting (Fig. 3G) illustrate that CellMentor preserves distinct clusters, whereas PCA collapses into mixed populations and Harmony shows only partial recovery.

Finally, we assessed the ability to handle rare cell populations by simulating a scenario where one of the cell types represented only 1% of cells across all batches (Fig. 3I). In this challenging setting, other tested dimensionality reduction methods failed to identify the rare population, with these cells appearing scattered throughout the embedding space with no distinct clustering. Standard approaches like PCA and even integration methods like Harmony were unable to preserve the signal from this rare population. CellMentor, however, leveraging its supervised framework and the presence of this rare population in the reference data, successfully maintained separation of this rare cell type in the query data, forming a distinct, coherent cluster.

Application to diverse tissue types
We next evaluated CellMentor on real datasets spanning multiple scenarios: paired reference-query datasets from specific tissues and comprehensive cross-technology transfer analysis. We first tested three pairs of datasets from different tissue types, each consisting of a reference dataset for training and an independent query dataset for validation: pancreatic tissue (Baron25/Muraro26), peripheral blood mononuclear cells (PBMCs27/Ding28), and melanoma samples (Tirosh29/Jerby-Arnon30). After standard preprocessing (Methods), we applied CellMentor and competing methods to each tissue type (Supplementary Fig. 9).

For pancreatic tissue analysis (Fig. 4A), we observed that PCA-based, scANVI, and CellMentor performed well in separating the major pancreatic cell types. This relatively high performance across these methods likely reflects the distinct transcriptional profiles of different pancreatic cell types, making this a more straightforward classification problem. During our analysis, we identified 16 cells that were annotated as alpha cells in the original Muraro dataset but were clustered with the beta cells according to CellMentor. We further investigated these cells using SingleR for independent annotation, which confirmed they expressed gene signatures more consistent with beta cells than alpha cells. This finding demonstrates how our approach can help identify potential misannotations in published datasets.

In PBMC analysis, the query dataset (Ding) contained samples from two different individuals, presenting a common challenge in single-cell analysis where batch effects often dominate technical variation (Fig. 4B). The PCA-based dimensionality reduction failed to integrate these batches, creating distinct clusters for each patient rather than grouping cells by biological type. While both Harmony and scANVI successfully integrated the batches, they showed limitations in preserving cell type structures—Harmony created multiple clusters for B cells and CD8+ T cells, while scANVI failed to clearly isolate natural killer (NK) cells. In contrast, CellMentor not only successfully integrated the batches but also generated well-defined clusters that corresponded more precisely to known cell types. Notably, none of the methods, including CellMentor, created completely distinct clusters for CD4+ T cells, CD8+ T cells, and NK cells. This likely reflects the biological reality that these cell types exist on a transcriptional continuum rather than as discrete populations, highlighting a limitation of discrete clustering approaches for certain cell types.

The melanoma datasets (Fig. 4C) presented the most challenging scenario, with complex mixtures of tumor and infiltrating immune cells that typically do not separate as cleanly as cells from normal tissues. This scenario represents many real-world tumor datasets where cellular states exist along continuums rather than as discrete types. PCA-based dimensionality reduction showed poor separation between similar cell types; CellMentor maintained clearer boundaries between different cell types while still effectively integrating cells from different patients.

Quantitative evaluation using the ARI across an expanded comparison with eleven other methods showed that CellMentor was among the top performers across all three tissue types (Fig. 4D and Supplementary Fig. 9A). CellMentor achieved an average ARI score of 0.710, virtually identical to other leading methods, including rssNMF (0.724), CASSL (0.723), and SIDA (0.710), and substantially outperforming traditional approaches like PCA (0.660) and other established methods. CellMentor also achieved a high average NMI score across the datasets (Supplementary Fig. 9B). These results demonstrate CellMentor’s consistent effectiveness across diverse tissue types and experimental settings, establishing it as a reliable method that performs at the state-of-the-art level alongside the best available alternatives.

To assess CellMentor’s ability to transfer knowledge across different sequencing platforms, we conducted a comprehensive analysis using datasets from Tabula Muris31, testing whether cell type-specific patterns learned from 10X Chromium reference datasets could effectively analyze Smart-seq2 query datasets across ten diverse tissue types: bladder, kidney, limb muscle, liver, lung, mammary gland, marrow, spleen, trachea, and tongue. CellMentor demonstrated superior cross-technology transfer performance, achieving the best performance among all methods in 5 out of 10 tissue types and remaining competitive in the remaining 5 tissues (Fig. 4E, F and Supplementary Fig. 6). Across all tissues, CellMentor achieved the highest average ARI score of 0.805, substantially outperforming the second-best method, PCA-based (Seurat), which achieved an average ARI of 0.671. This robust cross-platform performance demonstrates CellMentor’s ability to learn technology-invariant biological signatures that generalize effectively across major sequencing platforms. Performance variation across tissues likely reflects differences in cellular complexity and annotation quality, as these reference annotations were originally derived using traditional clustering approaches that may not capture the same cell type boundaries that supervised methods can identify.

Computational performance
Computational analysis revealed that CellMentor’s runtime is competitive with other NMF-based approaches across dataset sizes (Supplementary Fig. 7). While CellMentor can require longer runtimes when extensive hyperparameter tuning is performed, our systematic parameter sensitivity analysis identified robust default parameter ranges that provide strong performance without requiring exhaustive optimization (Supplementary Fig. 8). Using these evidence-based defaults significantly reduces computational requirements while maintaining CellMentor’s superior performance. The decomposition phase represents the most computationally intensive step, while the projection phase (applying to query data) is relatively fast. This makes the method practical for applications where the same reference will be used for multiple query datasets, as the computational investment in learning the reference representation is amortized across multiple applications. Our analyses show that memory usage scales linearly with dataset size, similar to other matrix factorization approaches. CPU time increases predictably with both cell number and the extent of hyperparameter optimization performed, but the availability of validated default parameters eliminates the need for extensive tuning in most applications.

Methods
CellMentor algorithm overview
CellMentor is a supervised dimensionality reduction framework designed specifically for scRNA-seq data that operates in two phases: decomposition (training) and projection (inference). The method uses labeled reference datasets to learn biologically meaningful latent spaces that can be effectively transferred to new datasets.

Decomposition phase: given a labeled reference dataset X ref 2 Rm×n with known cell type annotations, CellMentor employs supervised NMF to decompose the data into two matrices:

X ref (cid:2)W(cid:3)H ref (1)

where W 2 Rm×k represents gene combinations (meta-genes) and H ref 2 Rk×n contains cell-specific usage of these factors. The key innovation is that this decomposition is guided by cell type labels through discriminative constraints that explicitly preserve cell type identity while maximizing between-class separation.

The decomposition phase addresses three critical challenges:
1. Optimal rank selection: determining the number of factors (k) that capture biological variation while excluding technical noise through biwhitening analysis combined with eigenvector localization.
2. Matrix initialization: establishing robust starting points for the optimization.
3. Supervised learning: incorporating discriminative constraints that maximize between-class separation while minimizing within-class variation.

Projection phase: CellMentor uses the learned W matrix to analyze new query datasets. For a query dataset X query, we solve:

min H query ≥0 kW(cid:3)h query,j (cid:4) x query,j k2 (2)

This is solved using the Lawson–Hanson non-negative least squares algorithm. The resulting H query matrix represents the query cells in the biologically meaningful latent space learned from the reference.

Supervised non-negative matrix factorization. Since the core innovation of CellMentor lies in its supervised learning framework, we start with it. Our objective function balances five components:

F(W,H,α $, β $, γ, δ) = 1 kX(cid:4)WHk2 + S ~ W (cid:4) S ~ B + γ H + δ trace(WTWM) (3)

The reconstruction term 1 kX(cid:4)WHk2 ensures the factorization accurately represents the original data.

The within-class compactness term S ~ W minimizes variation among cells of the same type:

S ~ W = 1 trace((H(cid:4)HA)N ~ (H(cid:4)HA)T) (4)

where C is the number of cell types. This term computes how far each cell is from the average representation of its cell type in the latent space. The matrix A computes class averages, so (H−HA) represents the deviation of each cell from its class centroid.

The between-class separation term S ~ B maximizes distances between different cell types:

S ~ B = 1 trace((HBP)β $ (HBP)T) (5)

This term encourages different cell types to be well-separated in the latent space. The matrix B maps cells to their class identities, while P encodes all pairwise class differences, allowing us to compute distances between all cell type centroids simultaneously.

In addition, we include a sparsity constraint γ ∑ H i,j that promotes sparse factor representations, and an orthogonality constraint δ 2 trace(WTWM) that reduces redundancy between basis vectors, where M is a matrix with ones everywhere except the diagonal.

Cell-specific weights (α $ ): diagonal matrix where each element weights individual cells based on their distance from class centroids and class-level errors:

α $ (i,i) = d(cell i , centroid class(i)) × error rate class(i) (6)

Cells farther from their class center receive higher weights, encouraging the optimization to bring them closer to their class.

Class-pair weights (β $ ): diagonal matrix weighting class pair separations based on confusion frequency:

β $ (q,q) = confusion frequency(i, j) (7)

Class pairs that are frequently confused receive higher weights, increasing the penalty for insufficient separation. These weights are updated iteratively during optimization using intermediate SingleR predictions, allowing the algorithm to focus effort on the most challenging discrimination tasks.

Finally, the non-convex optimization problem is solved subject to non-negativity constraints W ≥ 0 and H ≥ 0 using multiplicative update rules:

W update:

W(t+1) = Wt (cid:5) WtHt(X)T + δWtM (8)

H update:

H(t+1) = Ht (cid:5) WTX + (WTX)2 + 4Φt (9)
2(WTWHt + Θt)

where Φt and Θt incorporate the discriminative constraint terms, and ⊙ denotes element-wise multiplication.

Optimal rank selection. Traditional rank selection methods fail for scRNA-seq data due to high noise levels and sparse structure. CellMentor uses a two-step approach:
Step 1: biwhitening analysis
We first test different scaling approaches (Poisson vs. quadratic variance) to normalize gene expression variance, then apply the biwhitening transformation to achieve unit variance in both dimensions. We compare the resulting empirical eigenvalue distribution to the theoretical MP distribution and select the scaling approach that minimizes the Kolmogorov-Smirnov distance to the MP distribution.
Step 2: eigenvector localization filter
Not all eigenvalues exceeding the MP threshold represent biological signal. We distinguish signal from noise using eigenvector localization, where localized eigenvectors have components concentrated in specific gene sets (indicating biological programs), while delocalized eigenvectors have Gaussian-distributed components (indicating random noise). For each eigenvector with an eigenvalue above the MP threshold, we test whether its components follow N(0, 1) using the Kolmogorov-Smirnov test. If we reject the null hypothesis (p < 0.01), the eigenvector is considered localized and represents a biological signal. If we accept the null hypothesis (p ≥ 0.01), the eigenvector is delocalized and represents technical noise.
Final rank: count of localized eigenvectors, effectively filtering technical noise while preserving biological dimensions.

Matrix initialization. Matrix initialization critically affects NMF convergence due to non-convexity. We implemented five strategies:
Random approaches:
(cid:129) Uniform: sample from [ϵ, max(X)] where ϵ is machine precision and max(X) is the maximum value in the dataset. Zero values are excluded from initialization to prevent elements from remaining permanently at zero during matrix updates.

Supervised approaches:
(cid:129) Up/down-regulated meta-gene (regulated): begin with uniform random initialization in the interval [ϵ, 1], then calculate the total sum of all H matrix elements and add this sum to the ith meta-gene for all cells belonging to the ith class (for i=1, 2, ..., min(k, C) where C is the number of classes). This creates class-specific emphasis by preferentially activating certain meta-genes for specific cell types, encouraging identification of up-regulated and down-regulated gene programs.

Structured approaches:
(cid:129) NNDSVD: non-negative double singular value decomposition that performs two sequential SVD operations to ensure non-negative initial matrices. The first SVD computes X = U Σ VT, while the second approximates the positive sections of the derived factors. Zero entries are replaced with avg(X) × Uniform([0, 1]) to avoid permanent zeros during optimization.
(cid:129) Gene clustering (skmeanGenes): initialize W using fuzzy c-means clustering to partition genes into k clusters, where each gene receives a degree of membership (between 0 and 1) to multiple clusters based on proximity to cluster centers. The membership degrees initialize the rows of W, and the H columns are initialized as the average projection values (WTX) for each cell type.
(cid:129) Cell clustering (skmeanCells): initialize H using k-means clustering to divide cells into k clusters, where the dominant cell type in each cluster determines the cluster identity. Cells not matching their cluster’s dominant type are reassigned to the most similar cluster of their actual cell type using cosine similarity. H columns are initialized with cell-cluster similarities, and W rows are initialized with cluster centroids.

After a comprehensive evaluation across multiple datasets, we found that the regulated meta-gene approach consistently produced the most robust results. This approach has been set as the default initialization method in CellMentor. However, all initialization strategies remain implemented in the package, allowing users to select alternative methods if desired for specific applications.

Hierarchical parameter optimization. CellMentor requires tuning four parameters (γ, δ, α, β). We use a hierarchical approach, separating global parameters:
Global parameters (γ, δ):
(cid:129) Grid search over: γ ∈ {0, 0.001, 0.01, 0.1, 0.5, 1}, δ ∈ {0, 0.1, 0.5, 1, 2, 5}.
(cid:129) Evaluate each combination using a train-validation split.
(cid:129) Select parameters maximizing validation clustering accuracy.

Notation and definitions. A summary of all symbols and notations used throughout the framework is provided in Table 1.

Simulation framework
We used Splatter (v1.26.0)17 to generate synthetic scRNA-seq datasets for systematic evaluation of CellMentor. Each simulation consisted of either two or four batches (1000 cells per batch) and six cell types with varying degrees of technical noise. Ten distinct simulation scenarios were created with progressively increasing difficulty to thoroughly test algorithm robustness. For each simulation, one or two batches were designated as reference data for the decomposition phase and the remaining batches as query data for the projection phase, mimicking real-world scenarios where reference and query data come from different experimental conditions.

The simulations were designed with systematically increasing difficulty (simulations 1–10) by manipulating key parameters. Early simulations (1–2) featured more balanced cell distributions, minimal batch effects (batch.facLoc=0.05, batch.facScale=0.1), low biological variation (bcv.common=0.2), and uniform dropout rates. The batch.facLoc and batch.facScale parameters control both the magnitude and fraction of genes affected by batch effects, directly testing scenarios where technical variation could confound cell type discrimination. Middle simulations (3–6) introduced stronger batch effects, higher biological variation (bcv.common increasing to 0.5), and heterogeneous dropout patterns across batches. The most challenging simulations (7–10) incorporated high biological variation (bcv.common=0.7–0.8), stronger batch effects (batch.facScale up to 0.5), higher outlier probability (out.prob=0.1), more heterogeneous dropout patterns, and highly skewed cell type proportions with up to 30% of cells in one type and as low as 3–7% in rare cell types. To assess statistical robustness, each simulation was repeated five times using different random seeds (seeds 100, 200, 300, 400, 500).

To test specific challenges, we designed additional simulation experiments exploring challenging scenarios in single-cell analysis. First, we evaluated novel cell type identification by excluding one cell type from the reference data while maintaining it in the query data. Second, we tested performance with incomplete reference annotation by merging two distinct cell types under a single label in the reference dataset, effectively creating corrupted reference data where annotation granularity was insufficient to capture biological subtypes. Third, we assessed robustness to reference annotation quality by systematically introducing mislabeling. For random mislabeling, we randomly reassigned 5%, 10%, 20%, 30%, 40%, and 50% of reference cells to incorrect labels without modifying their expression profiles. For correlated mislabeling, we simulated more biologically realistic annotation errors by modifying both labels and expression data. First, we computed centroid distances for all cells and identified boundary cells whose distance to the nearest incorrect cell type centroid was within a factor of 1.3 of the distance to their true centroid. From these candidates, we randomly selected cells to mislabel at rates of 10%, 30%, 50%, 70%, and 90%. For each mislabeled cell, we blended its expression profile with the centroid of its incorrect target type across up to 50 of the target’s top discriminatory genes, using a blending strength of 0.3. This process preserved realistic expression distributions while creating transcriptionally ambiguous cells that more closely resemble real-world misannotation scenarios where boundary cells are most susceptible to misclassification.

Fourth, we simulated technical variation affecting discriminatory features. We first identified the top 150 discriminatory genes for each cell type based on absolute fold change between the target population and all others. We then selectively perturbed a controlled fraction of these genes (25%, 50%, 75%, or 100%) across a controlled fraction of cell types (25%, 50%, 75%, or 100%). Perturbations were applied multiplicatively, introducing log-fold changes of approximately two (corresponding to 4-fold shifts in expression), with log-normal per-cell jitter (sdlog=0.25) to avoid perfectly coherent effects. Perturbations were introduced either in the reference batches only, the test batches only, or both simultaneously, allowing systematic evaluation of how corruption of the most informative features impacted downstream performance.

Fifth, we tested rare cell type detection by setting one cell type to represent only 1% of the total population across all batches.

All parameters used in the simulation framework are summarized in Supplementary Data 1.

Real data analysis
We evaluated CellMentor on three pairs of published scRNA-seq datasets spanning different tissue types, each consisting of a reference dataset for training and an independent query dataset for validation:
(cid:129) Pancreatic tissue: Baron human pancreas dataset (reference, n=8569 cells) and Muraro pancreas dataset (query, n=2042 cells). For feature selection, we supplemented variable genes with known pancreatic markers, including GCG, INS, SST, PPY, GHRL, KRT19, CPA1, and PDGFRB. In the Muraro dataset, mesenchymal cells were excluded to maintain consistency with other pancreatic datasets.
(cid:129) PBMCs: Zheng dataset (reference, n=6932 cells) and PBMC Single Cell Atlas dataset (query, n=7120 cells). Cell type-specific markers including CD19, CD79A, MS4A1, CD3D, CD8A, IL7R, and NKG7 were included in the feature list. In the PBMC dataset, we excluded cells with high mitochondrial gene expression (>5%) and abnormally high feature counts (>4000) as indicators of poor quality.
(cid:129) Melanoma: Tirosh melanoma dataset (reference, n=4097 cells) and Jerby-Arnon dataset (query, n = 6879 cells), with feature selection supplemented by known markers including CD2, CD3D, CD19, CD79A, and melanoma-specific genes such as MIA and TYR.

Data preprocessing included log-transformed normalization (ln(count × 10,000 + 1)), feature selection using the SingleR (v2.4.1) library; trainSingleR algorithm, and scaling using root-mean-square-deviation. For each tissue type, we retained only genes shared between reference and query datasets (pancreas: 15,290 genes; PBMC: 15,752 genes; melanoma: 22,667), selecting the most variable genes for analysis (pancreas: 755; PBMC: 487; melanoma: 562).

We furthermore evaluated performance across different sequencing technologies using the Tabula Muris consortium datasets, comparing 10x Chromium droplet-based sequencing with Smart-Seq2 plate-based sequencing across multiple mouse tissues. Tabula Muris data were obtained through the ExperimentHub package, comprising both droplet-based (10x Chromium) and Smart-Seq2 datasets. We identified tissues present in both collections and applied stringent quality control criteria: minimum 50 cells per dataset, at least 2 common cell types between technologies, and exclusion of cells with undefined cell ontology classifications. Additional quality filters removed cells with high mitochondrial gene expression (>5%) and abnormally high feature counts (>4000) as indicators of poor quality or doublets. For each tissue, we used the 10x Chromium dataset as a reference for training and the Smart-Seq2 dataset as a query for evaluation.

All analyses were performed on distinct single cells, with each cell representing an independent measurement. No repeated measurements of the same cells were conducted across any of the analyses.

Data analysis was performed using R (v4.3.2). CellMentor utilizes the following packages: methods (v4.3.2), Matrix (v1.6-5), sparsesvd (v0.2-2), RMTstat (v0.3.1), SingleR (v2.4.1), MLmetrics (v1.1.3), nnls (v1.5), parallel (v4.3.2), progress (v1.2.3), cluster (v2.1.4), skmeans (v0.2-17), irlba (v2.3.5.1), lsa (v0.73.3), data.table (v1.16.2), ggplot2 (v3.5.1), magrittr (v2.0.3), tibble (v3.2.1), Seurat (v5.1.0), utils (v4.3.2), matrixStats (v1.4.2).

Comparison with existing tools
We systematically compared CellMentor against thirteen dimensionality reduction methods spanning three categories: unsupervised approaches (PCA, ICA, NMF, pCMF, GLM-PCA), semi-supervised methods (CASSL, rssNMF, LIGER, scCoGAPS, scANVI, SiameseDA), integration-focused approaches (Harmony), and supervised discriminant analysis method (SIDA). All methods were implemented with consistent parameters (20 latent dimensions) and evaluated using standardized workflows. PCA was implemented using Seurat (v5.1.0)4 with 2,000 variable features. Independent Component Analysis (ICA) was performed using the fastICA package (v1.2-5.1)18 with tolerance set to 10−6 and the logcosh contrast function. NMF was executed using the NNLM package (v0.4.4) with Lee’s multiplicative update algorithm, 500 maximum iterations, and 10−6 convergence tolerance. pCMF20 was implemented with both zero-inflation and sparsity parameters enabled using the pCMF package (v1.2.1). Generalized Linear Model PCA (GLM-PCA) was applied directly to count data using the glmpca package (v0.2.0)21 with Poisson family modeling.

For semi-supervised methods, CASSL was implemented using default parameters with –p_missing 0.5 –max_dim 15. rssNMF was applied using the rssNMF package14 with sparsity parameter α=2 and graph regularization parameter β=2 and 100 maximum iterations. LIGER was executed via rliger, with variable gene selection thresholds var.thresh=2.0–2.5 and alpha.thresh=0.7–0.8 to ensure ≤2000 variable genes. iNMF was run with k=20, λ=100, and 10 iterations. Quantile normalization and Louvain clustering (resolution=0.2) were applied to aligned datasets. scCoGAPS (v3.99.3) was run on reference count matrices using 500 iterations and 20 latent patterns. Test cells were projected via projectR() using the amplitude matrix. For single-cell variational inference (scVI), we utilized the Python-based scvi-tools (v1.2.1)8 framework with three hidden layers (256, 128, and 64 nodes), batch correction, and 250 training epochs with a 10−3 learning rate. SiameseDA (SIamese Domain Adaptation) was run using a Python implementation adapted from the original pancreas benchmarking script from the paper24, wrapped for automated I/O. For each reference–test pair, we intersected common genes, applied library-size normalization and log transformation inside the SiameseDA pipeline, and reduced dimensionality with truncated SVD to 50 components (–pca_dims=50). A 1D-CNN siamese backbone was trained with a joint loss (categorical cross-entropy for a classification head and a contrastive loss for paired distances) with loss weight α=0.6, Adam optimizer (learning rate 1×10−5), batch size 512, and 10 epochs (–epochs=10). The learned embeddings (64-dimensional) for reference and test cells were exported and added to Seurat as a SIDA reduction. SIDA (v1.0) was adapted to single-cell data using the duplicate-view approach, where the normalized expression matrix was duplicated to create two identical views. Only genes common to both reference and query datasets were retained. Cross-validation (5 folds, random grid search with 6 grid points) optimized the sparsity tuning parameters (Tau) for each view. Models were trained with standardization enabled, weight=0.5 (balancing separation and association), convergence threshold=0.001, and a maximum of 20 iterations. Optimal discriminant vectors from reference data were used to project query cells via matrix multiplication, with final cell type assignment by nearest centroid using Euclidean distance.

CellMentor was applied using cross-platform transfer learning with automatic parameter optimization over α, β ∈ {1, 5, 10}, γ=1, and δ=0.1, selecting the best model.

For both simulated and real data, we assessed each method’s performance using a standardized workflow: (1) application of dimensionality reduction, (2) UMAP visualization of latent space, (3) graph-based clustering with multiple resolution parameters (0.1–0.5), and (4) evaluation of clustering quality against ground truth labels (simulations) or trusted annotations (real data). Clustering performance was quantified through ARI and NMI using the aricode package (v1.0.3).

Clustering and evaluation protocol
For all methods, clustering was performed on the full reduced-dimensional representation, not on 2D visualization coordinates. Specifically, we applied graph-based clustering using Seurat’s FindNeighbors and FindClusters functions on the latent space produced by each method. For PCA, ICA, NMF, pCMF, GLM-PCA, and scCoGAPS, we consistently used the first 20 dimensions (dims=1:20). For methods with data-driven or adaptive dimensionality—such as CellMentor, SCANVI, CASSL, and LIGER—the number of dimensions was extracted directly from model output (e.g., optimal factor number or latent bottleneck size), and the full latent representation was used for downstream analysis (typically 15–25 dimensions).

Clustering was performed using a fixed resolution parameter, chosen based on prior knowledge or method-specific characteristics (typically resolution=0.1 for PCA, ICA, NMF, and CellMentor; resolution=0.2 for methods such as SCANVI, LIGER, and scCoGAPS). No resolution grid search or optimization was conducted.

UMAP visualization was applied solely for interpretability, projecting the high-dimensional reduced space into two dimensions for display. All quantitative evaluations (ARI, NMI) compare cluster assignments from the full-dimensional clustering to ground truth cell type labels, ensuring that performance metrics reflect the quality of the reduced-dimensional representation rather than visualization artifacts.

Parameters sensitivity analysis. To assess CellMentor’s robustness to parameter selection, we conducted sensitivity analysis across the five key hyperparameters using simulated datasets with known ground truth (Supplementary Fig. 8). We tested α (within-class compactness) ∈ {0.1, 0.5, 1, 2, 5, 10, 20}, β (between-class separation) ∈ {0.1, 0.5, 1, 2, 5, 10, 20}, γ (sparsity) ∈ {0, 0.001, 0.01, 0.1, 0.5, 1}, δ (orthogonality) ∈ {0, 0.1, 0.5, 1, 2, 5}, and five initialization methods (uniform, regulated, NNDSVD, gene clustering, cell clustering), resulting in 12,005 parameter combinations. The analysis revealed that matrix initialization has the strongest influence on performance, with the regulated method consistently outperforming alternatives through its class-specific meta-gene emphasis. The sparsity (γ) and orthogonality (δ) parameters showed moderate sensitivity, while the discriminative parameters α and β demonstrated relatively low sensitivity despite their central role in the loss function. This low sensitivity of the discriminative parameters indicates that CellMentor’s supervised framework provides inherent stability once discriminative constraints are active, reducing dependence on precise parameter tuning.

Based on this systematic evaluation, we established evidence-based default parameter ranges: α ∈ {0.1, 0.5, 1}, β ∈ {1, 2, 5}, γ ∈ {0, 0.1}, δ ∈ {0, 0.5}, with regulated initialization as the recommended method. These defaults provide robust performance across diverse datasets without requiring extensive optimization for each application.

Integration with standard single-cell workflows
CellMentor is designed as a drop-in replacement for PCA in standard single-cell analysis pipelines. The typical workflow proceeds as follows:
Cluster-based annotation: apply standard clustering algorithms (e.g., graph-based clustering) to the CellMentor reduced dimensions, then annotate resulting clusters using established approaches such as marker gene analysis, cluster-level SingleR annotation, or manual curation based on known cell type signatures. CellMentor’s advantage lies in producing better-separated clusters that correspond more accurately to true cell types, making subsequent annotation more reliable.
Direct cell-level annotation: users can apply existing annotation methods like SingleR directly to the CellMentor latent space rather than to the original high-dimensional data, potentially benefiting from the denoised, biologically-focused representation.

The choice between approaches depends on specific analysis goals, with CellMentor’s strength being the creation of latent representations where cell types are maximally separable, benefiting both clustering and direct annotation workflows.

Computational performance analysis
We assessed computational performance across all methods using three dataset sizes: small (500 cells), medium (2000 cells), and large (5000 cells). All analyses were performed on identical hardware (Linux server running Ubuntu 20.04.6 LTS. The system was equipped with 4× Intel Xeon Gold 6248 CPUs (20 cores each, 160 threads total) and 1.97 TB of RAM. All methods were executed on the CPU. Analyses were conducted using R (version 4.3.2) within the RStudio environment or Python (version 3.12.2).) with five independent runs per method-dataset combination. We measured wall clock time, CPU time, and peak memory usage. For CellMentor, timing includes both the decomposition phase (learning from reference data) and the projection phase (applying to query data), with hyperparameter optimization using default parameter ranges.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "I thank the authors for kindly addressing my prior comments, and offering more clarity. I would still like a more thorough understanding of the benefits and limitations of the proposal, given the competition. I present my thoughts below.\n1. Data normalization: It is unclear to me if benchmarking for both simulations and real world datasets used the same preprocessed data. For example, page 16’s first paragraph says:”Data preprocessing included log-transformed normalization … feature selection using the SingleR (v2.4.1) [2] trainSingleR algorithm, and scaling using root-mean-square-deviation (RMSD). ”, and then in the last paragraph of the same page, I see “quantile normalization” being mentioned for iNMF. In the same paragraph, I see Liger operating on different gene selection thresholds. And so on. So I am confused : we need to benchmark the methods such that all the methods that allow the same preprocessed data input as CellMentor / options to generate the same data input as CellMentor, are setup such that they are all operating on the same data. Could the authors clarify this?\n2. Methods comparison: SIDA & DeepIDA : Considering that the central idea of these two methods is the same as CellMentor, the fact that they were developed for multi-modal data integration should not stop us from benchmarking it for single-cell data dimension reduction & cell type identification. Indeed, analysts “integrate” batches of single-cell data all the time (e.g., Harmony, Liger), and the multiple batches can be simply thought of as multiple modes of measurements. Furthermore, SIDA & DeepIDA methods reconstruct far fewer number of factors than the number of samples or features. Their low rank outputs can indeed be used for further cell type clustering. In the event that running SIDA & DeepIDA on large datasets is infeasible (potentially owing to the way they are currently implemented), comparisons with CellMentor should still be made under a fixed tractable number of features, with the computational considerations highlighted as point of impact. These comparisons will allow us to understand the impact of the contributions made in this paper.\n(Given my strong emphasis of these two techniques in my two reviews, I wish to declare in case it matters, that I was not involved in the development of SIDA/DeepIDA. I also do not know or collaborate with the group).\n3. I am still not convinced that we have understood the impact and limitations of the using discriminatory directions, in the presence of technical variation affecting discriminatory features, and limited biological knowledge in reference.\n3.1. As pointed out in my prior review’s comment 1: we still need to assess when cell type discriminatory features are perturbed by batch variation.\nSimulation of batch effects : Splatter unfortunately alters every gene by its own batch factor. This text from the manuscript “The batch.facLoc and batch.facScale parameters control both the magnitude and fraction of genes affected by batch effects, directly testing scenarios where technical variation could confound cell type discrimination. “ is unfortunately not true. These two parameters in Splatter instead reflect measures of mean and variance of the distribution from which the batch factors are estimated. There is nothing in Splatter that allows us to systematically induce batch effects in features (i.e., genes here) that discriminate cell types. Thus, there is no explicit understanding of what targeted perturbation the discriminatory features underwent in the various simulations. I therefore find that we still need to assess CellMentor’s performance as a fraction of cell type discriminatory features are progressively perturbed by batches.\nThis was suggestion 1.a of my comment in my previous review.\n3.2. Mislabeling experiment: Real world mislabeling of a cell occurs because its count data appears more closer to that usually seen from other types. In other words, merely swapping labels of generated count data is insufficient — the count distribution must also correlate with the labels. Because the current mislabeling experiment relabels portions of the simulated count data without correlated errors in the count data as well, I do not agree that the mislabeling experiment offers strong evidence of robustness. Decoupling count data correlated with misleading labels is a much harder problem to solve than what Supplementary Figure 5 would lead us to think.\nThis correlated patterns is what my comment 1.b from prior review was getting at, by simulating cell-types and subtypes of celltypes. The pancreatic example 4A with alpha /beta cells offer slight evidence to this point.\nIf we can show evidence for superior performance in 3.2 and 3.3, this would be a significant point of impact in favor of CellMentor.\n4. Minor Comments:\nA. Because we don’t know how a query dataset will behave, it is best to avoid the following sentence in page 5 (for example, think of the situation when both the reference and query are IID, in this case there won’t be any differences in how they project). As it stands, I find the claim here a bit of a stretch: “An interesting observation from comparing reference and query visualizations is that while reference datasets may retain some batch structure in the learned latent space, the projection phase effectively normalizes these effects when applying the learned patterns to query data (Supplementary Figure 4)”.", "reviewer_2": "All comments have been addressed and the manuscript substantially improved.\n(Remarks on code availability)\nI previously reviewed the code. It is usable, available, and ran with no issues.", "reviewer_3": "I appreciate the authors’ careful and thorough responses to my earlier feedback. I have just a few additional questions and points for clarification.\n1. Supp Fig 8A showed NMIs can be as low as 0, which is concerning. How many replicates were used? Are there anything special with these NMI runs — features in the input, intermediate results, or final outputs that may allow you to send a warning to users about the validity of the results?\n2. Regarding the matrix initialization method, it is not clear what “enhanced with class-specific meta-gene emphasis” means. Also, the terminology used in the methods section, Supp Fig 8B plot, and Supp Fig 8B legend are not consistent.\n3. According to Supp Fig 7 memory usage plot, “median size” of 2k cells needs more than 1e+05 MB (~100 GB) of memory. This is a very high requirement, especially considering the growing size of the datasets in recent publications. Also the PCA memory usage seems high. I’d recommend double check the statisitics or clarify how it was measured just in case.\n4. Other small issues (according to the PDF in the reviewer portal, maybe some of these are just formatting issues).\n- Supp Fig 8 was not included in the “SUPPLEMENTARY MATERIALS” section.\n- There are two section numbered as “4.1.1”.\n- Section “4.1.4 Notation and Definitions” is missing.\n- Section 4.7 was misplaced after the “Acknowledgement” section.\n(Remarks on code availability)\nI downloaded the code from https://github.com/petrenkokate/CellMentor on Aug 9th, 2025. I encountered exactly the same issues as I did during my first review when I run the code in the README file line by line (not skipping, adding, or switching any lines except executing \"library(scRNAseq)\" after I got the \"Error in counts(baron.h) : could not find function \"counts\"\" error message after running \"baron <- h.baron_dataset()\"."}, "author_response": "We thank the reviewer for raising this important point. All methods in our benchmark started from the same raw count matrices, and we applied a common preprocessing pipeline (library-size normalization, log1p, HVG selection, scaling) whenever the method could accept it. This ensured that PCA, ICA, NMF, and Harmony all ran on identical inputs with the same gene budget and identical train/test splits. Some methods, however, are defined on raw counts or require their own preprocessing. For example, CellMentor and scVI/SCANVI are designed to model counts directly, GLM-PCA is formulated for count input, pCMF requires its own filtering, and LIGER/iNMF has a specific normalization and gene selection procedure. In these cases, we followed the authors’ recommended pipelines but matched the gene budget across methods to keep comparisons fair. The mention of “quantile normalization” applies only to LIGER/iNMF and not to the other methods.\nFor full transparency, all preprocessing steps and parameters for each method are documented and reproducible in our benchmarking code, which is available on GitHub.\nWe thank the reviewer for highlighting this. We have added SIDA to our benchmarks under batch-effect simulations, and CellMentor consistently outperforms it (see updated Figure 2B and Supplementary Figure 3, which now also includes SIDA UMAPs).\nRegarding DeepIDA, we carefully considered whether it could be applied in our setting. However, DeepIDA was specifically designed for multi-omics bulk integration (e.g., transcriptomics, proteomics, metabolomics, lipidomics in COVID-19 cohorts). Its modeling framework assumes heterogeneous modalities with relatively small sample sizes, and its inference relies on correlating variation across these modalities. In contrast, single-cell RNA-seq data involves tens of thousands of cells from a single modality, with sparsity and dropout patterns that DeepIDA is not equipped to handle.\nEven under the reviewer’s suggestion of fixing a tractable number of features, the method remains conceptually mismatched: DeepIDA does not include cell-level modeling, clustering, or projection, but instead identifies cross-modal factors driving associations at the sample level. Thus, applying it to single-cell data would not yield a fair or interpretable comparison, as it is not designed to extract cell type–aware embeddings. For this reason, we believe benchmarking DeepIDA alongside CellMentor would be misleading, and instead focused our efforts on SIDA, which is applicable to single-cell scenarios.\nWe thank the reviewer for this critical observation about our previous simulation approach. The reviewer is absolutely correct that Splatter's batch parameters affect all genes globally rather than allowing targeted perturbation of discriminatory features. We have added a redesigned simulation framework to directly address this fundamental limitation.\nOur revised approach tackles exactly the scenario the reviewer described - what happens when technical variation specifically corrupts the genes that CellMentor relies on for classification. Instead of using Splatter's global batch effects, we developed a custom simulation that first identifies the 150 most discriminatory genes for each cell type using fold-change analysis, then selectively applies batch effects only to these critical features. This allows us to systematically test how performance degrades as an increasing fraction of discriminatory genes are corrupted by technical variation.\nWe implemented a comprehensive parameter grid that varies both the fraction of cell types affected by batch effects and the fraction of their discriminatory features that are perturbed, ranging from 25% to 100% in each dimension. This creates scenarios from minimal perturbation to maximum corruption of the discriminatory signals that methods depend on. We also tested three different scenarios for where these batch effects occur - in reference data only, test data only, or both simultaneously.\nThe results directly answer the reviewer's concern about CellMentor's robustness when discriminatory directions are compromised. Across all parameter combinations, CellMentor maintains substantially better performance than alternatives, achieving a mean ARI of 0.91 compared to 0.64 for PCA and 0.84 for Harmony. Most importantly, even in the most challenging scenario where all cell types have all their discriminatory features corrupted by batch effects, CellMentor maintains ARI above 0.8, while PCA drops below 0.2 and Harmony falls to around 0.4.\nWe added a new plot (Figure 3G-H) to demonstrate these results. The heatmaps in Figure 3H show CellMentor's consistent performance across the full parameter grid, contrasting sharply with the dramatic degradation of alternatives as discriminatory features are progressively corrupted. Figure 3G provides UMAP visualizations of the most challenging scenario, showing how PCA creates mixed clusters with poor separation, Harmony achieves partial recovery but with substantial cell type mixing, while CellMentor maintains clear boundaries between cell types even under maximum perturbation.\nThis targeted simulation framework provides the explicit understanding of discriminatory feature perturbation that the reviewer correctly identified as missing from our previous approach, directly demonstrating CellMentor's robustness when the features it relies on are systematically compromised.\nWe thank the reviewer for this insightful critique. We now added a new, redesigned, mislabeling simulation to address this fundamental issue. In our revised approach, we implement correlated mislabeling where cells that become mislabeled undergo actual transcriptional changes that make them genuinely resemble their new (incorrect) cell type. Specifically, we first identify boundary cells, by calculating distances from each cell to all cell type centroids and identifying cells that are naturally closer to other cell types than expected - these represent realistic candidates for annotation errors. Then, we create expression blending: For selected boundary cells, we blend their expression profiles with the centroid of their target (incorrect) cell type, focusing on the top discriminatory genes that define that cell type. This creates cells that express intermediate signatures between their true type and the mislabeled type. This blending process maintains realistic expression distributions while creating genuine transcriptional ambiguity that would fool both human annotators and computational methods.\nOur results show that this correlated mislabeling creates a much more challenging scenario, yet CellMentor was still able to perform reasonably. For example, at 50% mislabeling, CellMentor maintains an ARI of ~0.77 (Figure 3F).\nWe thank the reviewer for this comment. We understand that our statement was overly broad and could be misleading. We therefore revised this statement to be more cautious and acknowledge the limitations of our observation. The revised text reads: \"While our simulations suggest that the projection phase may help mitigate some batch effects present in the reference data (Supplementary Figure 4), this observation requires further validation across diverse experimental scenarios, as the effectiveness likely depends on the specific characteristics of both reference and query datasets.\"\nWe agree with the reviewer that NMI values close to zero are concerning and warrant clarification. These low values occurred in our sensitivity analysis runs, where certain simulation settings produced very sparse or uninformative features. We used 3 independent replicates for each parameter configuration. To improve transparency, we now added an automatic warning to the CellMentor package: if the resulting NMI falls below 0.2, the user is notified that the reference/test configuration may not contain sufficient discriminatory signal for meaningful annotation. This ensures that users are aware of potentially unreliable runs and can re-examine their input data or parameter settings.\nWe thank the reviewer for pointing out these issues with clarity and consistency. Regarding the unclear terminology, we agree that \"enhanced with class-specific meta-gene emphasis\" was vague and revised this to provide a clear description: the \"regulated\" initialization method begins with uniform random initialization in the interval [ε, 1], then calculates the total sum of all H matrix elements and adds this sum to the ith metagene for all cells belonging to the ith class (for i = 1, 2, ..., min(k, number of classes)), creating class-specific emphasis by preferentially activating certain metagenes for specific cell types. Following the comment, we recognized that our explanations of the matrix initialization strategies were incomplete and added more detailed descriptions for all five initialization methods to improve clarity and completeness.\nWe also verified the supplementary figure references to ensure they match the current manuscript structure.\nWe thank the reviewer for pointing out this discrepancy. Upon re-examining our benchmarking pipeline, we identified that the memory profiling function occasionally captured transient allocations from R’s garbage collector and intermediate objects, leading to inflated values in the original plot. We have corrected the profiling to consistently report peak resident memory used by each method, and we regenerated Supplementary Fig. 7 with the corrected values. In the updated figure, the memory usage for 2k cells is within expected. Also, we added SIDA method to our pipeline and to computational cost analysis accordongly.\nFixed.\nWe thank the reviewer for carefully checking the reproducibility of the README code. The error was caused by reliance on the counts() accessor without explicitly ensuring that the appropriate Bioconductor packages were loaded. In the revised version of CellMentor we modified the implementation of h.baron_dataset() to use fully qualified function calls from scRNAseq and SummarizedExperiment, and we ensured that these packages are listed under Imports in the DESCRIPTION file. As a result, the function now runs correctly without requiring the user to manually attach scRNAseq. To make the workflow more robust, we also added a note in the README recommending that users load the relevant Bioconductor packages manually if they encounter issues, although this should no longer be necessary in a clean installation."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I thank the authors for kindly addressing the reviewers’ comments. The effects of noise induced by mislabeling of cells in the reference, and batch variation on discriminatory features are now more reasonably assessed. CellMentor succeeds in effectively capturing cell type information.\nBelow I argue why, despite the authors' arguments, I still do not find SIDA & DeepIDA as conceptually mismatched tools. I ask again that we compare SIDA & DeepIDA. If running DeepIDA is complicated, we need to at least benchmark against SIDA.\nParenthetically, I must note that the names of tools have created confusion. I was referring to SIDA : Safo, Min and Haine, 2021, Biometrics, not the Siamese Domain Adaptation used by the authors. I think I had given the citation explicitly in the first review, but not emphasized it again in my second review. My apologies.\nFor reasons listed below, I find that SIDA & DeepIDA are not entirely conceptually mismatched.\n1. Both SIDA & DeepIDA are attempting to jointly classify groups of samples (e.g., patients), and find associations among features measured in the distinct modes of measurements (e.g., genes via bulk RNASeq, microbes via 16S rRNASeq say). The output is a set of discriminatory vectors distinguishing groups (can be used as a low rank basis for samples, our main interest), and across-mode association vectors for features (we don't care about this for this work).\n2. I realize both the tools are built for multimodal data integration. But any multimodal data integration can be run for single data analysis as well, by simply duplicating the single dataset into two copies and pretending that these are two modes. The fundamental goal here is to (A) obtain the basis for samples (let’s call it U, an n x K dimension matrix output by SIDA where n is the total number of samples and K is the total number of components), (B) so that we can obtain the feature loadings vector V (p x K, where p is the total number of features), and (C) use V to project a given query dataset. Because every measured feature is approximated as a linear combination of U, the corresponding feature loadings vector V is simply calculated by regressing every feature against U (example, see sparse PCA paper ). Thus, one obtains a pair (U,V) for every dataset from SIDA, from a reference dataset. Following the exact same logic, given a query dataset, one can project the query data using V via least squares $( V’V )^(-1) V’ QueryData$.\nThese steps are similar to several regularized PCA type approaches, and not new.\n3.Thus, if we pretend the known cell types in the reference as “groups”, and run the strategy in point 2 outlined above, we will be able to use SIDA to identify discriminatory directions, construct feature loadings vector V via regression, and project query data onto the components space via regression on V again. The authors could then assess the performance of SIDA’s output for the various benchmarks.\nMinor comment\n1. Where SIDA & DeepIDA will greatly suffer relative to CellMentor is in their additional goal of identifying across-modal association vectors of features (p-dimensional feature set), which necessitates working with large across-modal feature x feature covariance matrices. We shouldn’t have to do that if our focus is only single dataset analysis, and this is clearly an advantage for CellMentor.", "reviewer_3": "The authors have addressed my previous concerns in the revised manuscript. I am satisfied with the changes and have no further comments."}, "author_response": "We sincerely thank the reviewer for clarifying their reference to SIDA and for providing the detailed methodological framework for adapting it to single-mode data analysis.\nFollowing the reviewer's explanation, we implemented and extensively tested SIDA using the duplicate-view approach exactly as outlined.\nWe implemented SIDA following the three-step procedure the reviewer specified: duplicating each dataset to create two identical \"views,\" running SIDA on duplicated reference data with cell types as groups, calculating feature loadings through regression, and projecting query data. We used the official SIDA R package with default parameters and applied the same evaluation framework used for all other methods.\nOur evaluation revealed fundamental problems with SIDA's adaptation to single-cell analysis.\nWhen we ran SIDA 5 times with different random seeds for each of our 10 simulation scenarios (50 runs total), performance varied wildly from ARI 0.15 to ARI 1.00 across runs. This represents a standard deviation of 0.34 across all simulations, contrasting sharply with CellMentor's SD below 0.05 on the same data. The instability appears to stem from optimization challenges when duplicated views create perfect correlation, leading the alternating optimization to converge inconsistently to different local minima.\nThis variability creates severe practical problems. Results cannot be trusted even when SIDA appears successful, since we cannot distinguish genuine signal capture from favorable random initialization. Users would need multiple runs and manual result selection, which is methodologically problematic. The unpredictability makes SIDA unreliable for routine single-cell analysis.\nReal dataset testing confirmed these concerns. On the three initial tissue datasets (Pancreas, PBMC, Melanoma from Figure 4A-D), SIDA showed comparable performance to CellMentor, though with 16-28× longer runtime and doubled memory usage. However, the cross-technology transfer analysis across ten diverse tissues (Figure 4F) revealed critical reliability issues. Most concerning was performance on datasets like Trachea tissue, where SIDA achieved only ARI 0.37 while most other methods exceeded 0.6 and CellMentor achieved 0.77 (the highest among all methods). This catastrophic failure occurred on a relatively straightforward dataset, demonstrating unpredictable failure modes. This pattern repeated across multiple tissues, with SIDA failing to achieve competitive performance in 4 out of 10 tissues while CellMentor achieved the highest average ARI of 0.805 across all tissues. This unpredictability prevents us from predicting which datasets trigger these failures, rendering SIDA unreliable for general use.\nThe reviewer correctly identified that SIDA would suffer from computing large across-modal covariance matrices. Our results confirm this and reveal additional fundamental issues. The perfect correlation between duplicated views creates degenerate optimization landscapes. Single-cell datasets with tens of thousands of cells and genes require computing twenty-thousand by twenty-thousand covariance matrices twice, becoming computationally prohibitive. The duplication destroys natural sparsity in scRNA-seq data (typically 90% zeros), preventing efficient sparse matrix operations. Most fundamentally, SIDA was designed to find associations between different measurement types, but when both views are identical, the across-modal association component becomes meaningless and the optimization loses biological grounding.\nThe contrast highlights why purpose-built methods outperform adapted approaches. CellMentor optimizes directly on single-mode data without artificial duplication, maintains natural sparsity for efficient operations, exhibits stable optimization with minimal seed sensitivity, scales memory linearly rather than quadratically, and incorporates biwhitening-based rank selection tailored to scRNA-seq noise characteristics. SIDA's duplicate-view adaptation instead requires data duplication doubling memory, shows unstable optimization, forces dense matrix operations destroying sparsity, provides no scRNA-seq-specific noise handling, and applies a multi-modal framework to single-mode problems.\nWe have added SIDA comparisons throughout the manuscript: Figure 2 now includes SIDA in both simulation scenarios (with and without batch effects), Figure 4D and 4F include SIDA performance on real datasets, and we added comprehensive SIDA visualizations to Supplementary Figures 2, 3, and 9 showing both performance metrics and UMAP embeddings across all tested scenarios. We must first clarify an important naming issue to avoid confusion. In our previous response, we referenced \"SIDA\" when discussing Siamese Domain Adaptation methods. However, the reviewer correctly pointed out they were referring to SIDA from Safo, Min and Haine (2021, Biometrics) - Sparse Integrative Discriminant Analysis. To prevent further confusion, we have renamed the Siamese Domain Adaptation method to \"SiameseDA\" throughout our manuscript. All comparisons and discussions below refer exclusively to SIDA (Safo et al., 2021, Biometrics) as the reviewer intended.\nOur evaluation demonstrates that multi-view integration methods adapted through data duplication suffer from fundamental limitations for single-cell analysis. SIDA's extreme seed sensitivity (ARI 0.15-1.00 on identical data), catastrophic failures on standard datasets, and 16-28× longer runtime create insurmountable practical barriers. These findings validate the reviewer's observation that SIDA would greatly suffer but reveal issues more severe than anticipated.\nCellMentor's purpose-built design addresses these limitations through direct optimization on single-mode structure, native sparse operations, stable performance, scRNA-seq-specific rank selection, and atlas-scale computational efficiency. We have added SIDA comparison results throughout the manuscript and revised the text accordingly. This strengthens our work by demonstrating that successful supervised dimensionality reduction for single-cell genomics requires designs directly addressing scRNA-seq characteristics rather than adaptations of multi-modal frameworks."}, {"round": 3, "reviewer_comments": {"reviewer_1": "I thank the authors for kindly providing my requested comparisons to SIDA. As summarized succinctly by the authors, it is clear that besides the simple computational bottleneck, there are a variety of subtle and complicated issues in adapting SIDA developed for multi-modal data integration to single dataset analysis in single cells.\nI congratulate the authors on this beautiful work, and thank them for working with me through this process!"}, "author_response": ""}]

Decision letter