
==== Front
Proc Natl Acad Sci U S A
Proc Natl Acad Sci U S A
PNAS
Proceedings of the National Academy of Sciences of the United States of America
0027-8424
1091-6490
National Academy of Sciences

39226348
202400002
10.1073/pnas.2400002121
research-articleResearch ArticlestatStatistics433
Physical Sciences
Statistics
Single-cell analysis via manifold fitting: A framework for RNA clustering and beyond
Yao Zhigang zhigang.yao@nus.edu.sg
a 1 2
Li Bingjie a 1
Lu Yukun a 1 https://orcid.org/0009-0005-8617-3045

Yau Shing-Tung styau@tsinghua.edu.cn
b 1 2 https://orcid.org/0000-0003-3394-2187

aDepartment of Statistics and Data Science, National University of Singapore, Singapore 117546, Republic of Singapore
bYau Mathematical Sciences Center, Jingzhai, Tsinghua University, Beijing 100084, China
2To whom correspondence may be addressed. Email: zhigang.yao@nus.edu.sg or styau@tsinghua.edu.cn.
Contributed by Shing-Tung Yau; received January 29, 2024; accepted July 19, 2024; reviewed by Stephan Huckemann and Jingyi Jessica Li

1Z.Y., B.L., Y.L., and S.-T.Y. contributed equally to this work.

3 9 2024
10 9 2024
3 9 2024
121 37 e240000212129 1 2024
19 7 2024
Copyright © 2024 the Author(s). Published by PNAS.
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This open access article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

Significance

Single-cell RNA sequencing (scRNA-seq) analysis, crucial for uncovering cellular diversity and disease mechanisms, faces challenges due to technical variability, data complexity, and biological noise. Current scRNA analysis approaches often fail to deliver accurate and stable scRNA-seq clustering. In response, we present a pioneering framework centered on manifold fitting. Diverging from traditional methods that aim for a low-dimensional data representation, our approach fits a low-dimensional manifold within the ambient space and unfolds the data accordingly. This strategy effectively reduces distances between similar cell types while preserving expression information, significantly improving scRNA-seq clustering and visualization compared to existing state-of-the-art techniques. This groundbreaking advancement establishes a standard in scRNA-seq analysis, opening exciting avenues for future research and potential clinical applications.

Single-cell RNA sequencing (scRNA-seq) data, susceptible to noise arising from biological variability and technical errors, can distort gene expression analysis and impact cell similarity assessments, particularly in heterogeneous populations. Current methods, including deep learning approaches, often struggle to accurately characterize cell relationships due to this inherent noise. To address these challenges, we introduce scAMF (Single-cell Analysis via Manifold Fitting), a framework designed to enhance clustering accuracy and data visualization in scRNA-seq studies. At the heart of scAMF lies the manifold fitting module, which effectively denoises scRNA-seq data by unfolding their distribution in the ambient space. This unfolding aligns the gene expression vector of each cell more closely with its underlying structure, bringing it spatially closer to other cells of the same cell type. To comprehensively assess the impact of scAMF, we compile a collection of 25 publicly available scRNA-seq datasets spanning various sequencing platforms, species, and organ types, forming an extensive RNA data bank. In our comparative studies, benchmarking scAMF against existing scRNA-seq analysis algorithms in this data bank, we consistently observe that scAMF outperforms in terms of clustering efficiency and data visualization clarity. Further experimental analysis reveals that this enhanced performance stems from scAMF’s ability to improve the spatial distribution of the data and capture class-consistent neighborhoods. These findings underscore the promising application potential of manifold fitting as a tool in scRNA-seq analysis, signaling a significant enhancement in the precision and reliability of data interpretation in this critical field of study.

manifold fitting
single-cell RNA sequencing analysis
unsupervised clustering
visualization
Singapore MOE Tier 2 A-0008520-00-00 Zhigang Yao Singapore MOE Tier 2 A-8001562-00-00 Zhigang Yao Singapore MOE Tier 1 A-0004809-00-00 Zhigang Yao Singapore MOE Tier 1 A-8000987-00-00 Zhigang Yao
==== Body
pmcSingle-cell RNA sequencing (scRNA-seq) (1, 2) has become a crucial tool in genomic research, offering unparalleled resolution in dissecting the genomic, transcriptomic, and epigenomic profiles of individual cells. This detailed view is instrumental in deciphering the complex interplay within tissues and the inherent diversity of cellular populations. Advanced analytical methodologies (3, 4), encompassing dimensionality reduction, cell clustering, and visualization techniques, enable a deeper understanding of cellular growth tracks and the variations in gene expression. scRNA-seq has provided insights into the pathogenesis of diseases such as diabetes (5), Alzheimer’s disease (6), and cancer (7). Additionally, the advancement of scRNA-seq technology has set the basis for the growth of multiomics analysis (8), spatial transcriptomics (9), and the Human Cell Atlas project (10). These advancements not only enhance the depth of cellular analysis but also place these findings in the broader context of tissue structure and function.

Despite the significant contributions of scRNA-seq in genomic research, it faces notable challenges, particularly in dealing with two types of noise: biological noise and measurement error. Biological noise is intrinsic to the cells and can stem from various sources. This includes the overexpression of certain genes, nutrient fluctuations, the specific location of a cell within its tissue or organ environment, and the current state or condition of the cell. In contrast, measurement errors are tied to the technical aspects of the sequencing process. These errors can vary based on the sequencing technology and platform used, as well as the specific methods applied during sequencing. The combination of biological noise and measurement errors adds a layer of complexity and variability to scRNA-seq data, making it challenging to extract accurate and meaningful biological insights.

Researchers have developed various strategies to mitigate the effects of noise and variability in scRNA-seq data. These strategies can be broadly classified into three categories: genomic imputation, graph-based methods, and deep learning networks. Genomic imputation primarily addresses dropout events in scRNA-seq data, where a gene is expressed but fails to be detected. A notable example is CIDR (11), which calculates the probability of such dropout events and uses this information for imputation. Other genomic imputation methods, including MAGIC (12), SAVER (13), and scImpute (14), have been developed to address dropout events and other forms of noise, each employing unique strategies like Markov process, Bayesian prediction, and machine learning techniques. Graph-based methods, such as Seurat (15), utilize principal component analysis followed by graph clustering to categorize cells. Similarly, SC3 (16) employs a graph-based approach but focuses on consensus clustering. SCANPY (4) constructs a neighborhood graph of cells based on their similarity and then applies Louvain clustering (17) to this graph. Deep learning algorithms, like scDHA (18) and DESC (19), leverage computational models to represent and analyze scRNA-seq data. scDHA combines a nonnegative kernel autoencoder and a Bayesian autoencoder for clustering. DESC uses a stacked autoencoder for data representation and integrates it with an iterative clustering neural network.

While the strategies mentioned above provide some solutions to the noise issue in scRNA-seq analysis, each has limitations that can impact their effectiveness. Genomic imputation methods rely on assumptions about data distribution and dropout events, introducing the risk of biases or inaccuracies. Graph-based clustering methods, which involve dimensionality reduction, may lead to information loss, potentially obscuring key cellular differences crucial for understanding biological processes. Although deep learning algorithms are often more accurate, their decision-making processes may lack clarity. Furthermore, the specificity of deep networks could limit their adaptability to varied datasets.

Manifold fitting (20–23, 25, 26) is an advanced technique with the potential for processing the scRNA-seq data. This method aims to reconstruct a smooth manifold within the original space where the data are measured, capturing the low-dimensional structure of the data in a manner that minimizes information loss and effectively eliminates noise. The latest innovations in manifold fitting (23) offer solutions to the limitations of existing methods, characterized by three key features. First, it operates in the ambient space without the information loss typically associated with dimensionality reduction in graph clustering. Second, manifold fitting is adaptable to diverse data distributions, employing flexible neighborhood definitions, which is a distinct advantage over interpolation methods that rely heavily on specific data distribution assumptions. Finally, this technique is highly interpretable and supported by comprehensive theoretical analysis.

Building upon these advancements, we invent a framework named Single-cell Analysis via Manifold Fitting (scAMF) to address the persistent noise issue for scRNA-seq data. scAMF, inspired by recent advances in ref. 23, is elegantly designed to accommodate the analysis of scRNA-seq data. Compared with existing scRNA-seq analysis frameworks, scAMF carries several notable advantages. First, scAMF markedly enhances the spatial distribution of data, achieving a more pronounced agglomeration within classes while ensuring clearer separation between them. This refinement is pivotal for enabling more precise and accurate clustering in subsequent analyses. Second, scAMF comprehensively integrates various data transformation methods and clustering algorithms. This synergy allows for the efficient processing of clustering and visualization tasks across different scRNA-seq data platforms, demonstrating superior performance compared to existing algorithms. Most notably, scAMF introduces an innovative self-supervised approach to autonomously determine optimal clustering outcomes, marking a significant stride toward auto–machine learning in scRNA-seq analysis. Considering these advancements, scAMF emerges as a promising and potentially transformative tool in scRNA-seq analysis, setting a benchmark for future research.

Before diving into the details, we showcase the performance of scAMF in clustering on the Kolodziejczyk data. This set of mouse embryo stem cells (27) contains 704 cells in 3 classes (lif, 2i, and a2i) with 38,658 genes. This study focuses on the potential molecular mechanisms governing stem cell differentiation and maintenance. The true labels here in the Kolodziejczyk data represent the three culture conditions of mouse embryonic stem cells. These classifications are validated by biologists using gene markers, cell morphology, and cell location, similar to all other labels in our selected datasets used in the following sections. Fig. 1 shows the t-distributed stochastic neighbor embedding (T-SNE) visualizations (28) between raw and scAMF-processed data. While the data points belonging to different classes in the raw data are entangled together, scAMF manages to separate classes and form a distinct aggregated pattern within each class. This evidently suggests a much higher clustering accuracy, highlighting the great potential of scAMF for scRNA-seq analysis.

Fig. 1. An intuitive illustration of scAMF’s performance on the Kolodziejczyk data using T-SNE. The data points before and after manifold fitting are both reduced to two dimensions and colored according to the true cell types. After manifold fitting, the data points are visibly aggregated into three distinct clusters that align with the true cell types, rather than being entangled together without discernible clusters as observed before the fitting process. The clustering accuracy of k-means on the raw data is only 57%. However, after applying scAMF to fit the manifold, the accuracy increases to perfect 100%. The clear separation and the increase in accuracy evidently highlight the great potential of scAMF for scRNA-seq analysis.

In the following sections, we will explore scAMF further in several aspects on 25 scRNA-seq datasets (Results), with detailed procedure and methodology (Materials and Methods).

Results

An Overview of scAMF.

scAMF consists of three modules for the analysis of scRNA-seq data: raw data transformation, manifold fitting, unsupervised clustering and validation (Fig. 2A). Integrating these modules makes scAMF a unified framework for automatic machine learning of scRNA-seq data with high accuracy.

Fig. 2. The scAMF framework for the analysis of scRNA-seq data. (A) A schematic overview of the scAMF pipeline: raw data transformation, manifold fitting, and unsupervised clustering and validation. (B) An illustration of the mechanism of value-to-rank, unit-vector, and logarithmic transformations. The colored bars demonstrate randomly generated example sequences and the corresponding transformed results. The scatter plots depict the mean-variance distribution of the raw and transformed Kolodziejczyk data (27). (C) Illustration of neighborhood selection based on the shared nearest neighborhood: Each cell is denoted as ci, with its k-nearest neighborhood represented by Nk(ci). The shared nearest neighborhood pij quantifies the overlap between Nk(ci) and Nk(cj), indicating common points. This method refines the neighborhood of ci by selecting cells with the highest shared nearest neighborhood count. This process results in a more accurate neighborhood definition for ci.

scAMF first transforms the data using three transformation methods (Fig. 2B and Raw data transformation). The role of this module is to improve the signal-to-noise ratio in the data, including reducing excessive variance in highly expressed genes and partially correcting for batch effects, ensuring a more accurate and reliable data representation. Due to the diverse origins and sequencing platforms of scRNA-seq, different transformations are required to process the data effectively. Note that scAMF employs an unsupervised approach to determine the most suitable transformation method in its pipeline. This strategy is innovative among existing scRNA-seq analysis methods.

For each transformed scRNA-seq data, scAMF performs the manifold fitting algorithm [Manifold fitting, (23)], with suitable modifications tailored to scRNA data. To highlight, scAMF uses neighborhoods determined by the shared nearest neighbor metric (29), which is more efficient in the measurement of high-dimensional data than the classical Euclidean or correlation-based metrics. A demonstration of the principle that uses the shared nearest neighbor algorithm to determine neighborhoods is shown (Fig. 2C). Moreover, scAMF simplifies the two steps outlined in the original approach by Yao et al. (23), rendering it a more straightforward and expeditious method. Following the fitting process, the data are unfolded in the ambient space. This results in a reduction of intraclass spacing and an augmentation of interclass spacing, creating conditions more conducive to effective clustering in subsequent analyses.

For the fitted data, scAMF employs four fundamental and fast clustering algorithms, which include but are not limited to spectral clustering and three agglomerative hierarchical clustering algorithms with different linkage settings to accommodate various structural complexities within scRNA-seq datasets (Unsupervised clustering and validation). The utilization of multimodel clustering in scAMF is connected to SC3 (16). By combining these clustering techniques, scAMF ensures a comprehensive analysis, enabling the identification of subtle and distinct clustering patterns that a single-method approach might overlook. However, a key distinction from SC3 is that we do not rely on the consensus of multiple clustering results. Instead, we select the optimal clustering result, as introduced subsequently.

The unsupervised clustering validation module (Unsupervised clustering and validation) aims to identify the most accurate clustering result among outcomes produced by various transformations and clustering methods. This module operates on the principle that cells of the same type exhibit greater similarities while heterogeneous samples do not. For each clustering result, scAMF calculates the similarity ratio of intraclass cells to interclass cells. The clustering result with the highest similarity ratio is then selected as the most accurate one. Consequently, this process also determines the optimal transformation and manifold-fitting representation of the data.

By implementing scAMF, we can classify high-variance genes into two distinct groups: cluster-irrelevant genes and cluster-related genes. Cluster-irrelevant genes showcase variations not linked to specific cell types, essentially representing “housekeeping” genes essential for basic cellular functions but not indicative of particular cell types. In contrast, cluster-related genes exhibit significant variance across different classes, reflecting genes more indicative of specialized cellular functions and characteristics. scAMF identifies cluster-related genes based on clustering results, and these genes can be visualized. This facilitates a more focused and informative visualization, highlighting genes most relevant to the specific cell types and functions identified in the scRNA-seq data.

scAMF Enhances the Performance of scRNA-seq Clustering.

To evaluate the performance of scAMF in clustering, we compare it to seven cutting-edge scRNA analysis methods, including scDHA (18), SC3 (16), SINCERA (30), CIDR (11), Seurat (31), DESC (19), and SCANPY (4), on 25 scRNA-seq datasets with known cell types (Materials). Note that the true cell type is solely used for performance evaluation.

In evaluating clustering performance, scAMF demonstrates a notable superiority over other methods. We primarily employ the adjusted rand index (ARI) as our main evaluation metric, complemented by the normalized mutual information (NMI) and accuracy (ACC) metrics (SI Appendix, B). Here, ARI, NMI, and ACC are calculated between the predicted clustering labels and the ground truth labels of scRNA-seq data. The higher these indices are, the better the clustering algorithm has identified the inherent grouping. scAMF consistently registers the highest ARI scores across a majority of the datasets, averaging 0.78. This performance significantly eclipses that of its closest competitor, scDHA, which achieves an average score of 0.61 (Fig. 3A and SI Appendix, Table S2). The statistical significance of the superior ARI of scAMF is further validated by a one-sided Wilcoxon signed-rank test, which yields a P-value of 6.665×10−4. Additionally, scAMF’s exemplary performance in clustering is similarly evident in the NMI and ACC metrics (SI Appendix, Tables S3 and S4).

Fig. 3. Performance of scAMF on 25 scRNA-seq datasets. (A) Clustering performance comparison of scAMF and other methods, measured by ARI. The 25 plots display ARI values for each dataset, while the Lower Right plot shows the average ARI, with scAMF surpassing other clustering methods. (B) Scatter plots of ARI against UVI for all combinations of transformation methods and clustering algorithms in each dataset with manifold fitting. The blue circle represents the point with the highest UVI (hence selected by the scAMF), while the red circle represents the point with the highest ARI calculated by the true labels. In most cases, combinations with the highest UVI values align with those exhibiting the highest ARI, underscoring the rationality of scAMF’s selection procedure. (C) ARI of scAMF with and without manifold fitting. The inclusion of manifold fitting consistently enhances scAMF’s clustering performance.

The strong performance of scAMF is thanks in part to its unsupervised clustering validation module. It selects the clustering result that maximizes the unsupervised validation index (UVI) (Unsupervised clustering and validation) from combined transformation methods and clustering algorithms. Notably, the chosen result aligns with the highest clustering accuracy. This relationship is evident by comparing UVI and ARI across all clustering results for each dataset (Fig. 3B). In 16 out of 25 datasets, the highest UVI coincides with the highest ARI. Additionally, in 6 datasets, the clustering indicated by the peak UVI achieves the second-highest ARI. This highlights scAMF’s ability to effectively identify the optimal clustering result, contributing to automated machine learning with high accuracy. In instances where the ARI and UVI metrics do not align, we have carried out an in-depth examination to uncover the underlying factors. More information is available in SI Appendix, G and Figs. S1 and S7.

We conduct ablation studies to highlight the critical role of manifold fitting in scAMF. The impact of this step is illustrated by comparing ARI for scAMF with and without manifold fitting, as depicted in Fig. 3C. The results reveal a significant decrease in ARI; specifically, the ARI of scAMF without manifold fitting model drops by an average of 0.28 across nearly all datasets tested. This decline highlights the pivotal role of manifold fitting in enhancing scAMF’s performance. Additionally, the effectiveness of the unsupervised clustering validation module is considerably reduced without manifold fitting (SI Appendix, Fig. S2). In scenarios where manifold fitting is excluded while retaining other steps, only 7 datasets show a consistent correlation between UVI and ARI. This evidence collectively underscores manifold fitting as a cornerstone module in the scAMF framework, integral to its effectiveness.

scAMF Offers a Spatial Layout That Better Supports Class Separation.

To intrinsically evaluate the improvement of scAMF brought to scRNA-seq datasets, we employ specific metrics to evaluate spatial distributions post-scAMF application, focusing on intra/interclass distances and neighborhood purity.

Intraclass distance calculates the average distance from a sample to others within its class, indicating the compactness of each class. Conversely, interclass distance measures the average distance to a sample from different classes, illustrating class separation. Scatter plots across six typical datasets (Fig. 4A) show a trend toward increased class distinction under scAMF, indicated by the upper-left movement of data points. It visually confirms scAMF’s role in improving spatial class separation.

Fig. 4. Space distribution comparison between raw data and manifold-fitted data. (A) Representations of intraclass distance against interclass distance plots in six datasets: Kolodziejczyk, Goolam, Yan, Usoskin, Zeisel, and Darmanis. Fitted datasets exhibit higher intra–inter-class distance ratios, indicating improved underlying structures for clustering post-scAMF. (B) Box charts of raw and fitted data neighborhood purity for all 25 datasets. The majority of fitted datasets show higher neighborhood purity values than raw datasets, supporting the notion that scAMF contributes to a more accurate representation of scRNA-seq data structure.

The neighborhood purity of a sample is defined as the proportion of neighboring samples of the same type as the sample itself. It can be assessed before and after the scAMF framework. In our analysis, we choose a neighborhood size of ten and compute the neighborhood purity for all datasets. By comparing the results, we observed a consistent increase in neighborhood purity in the datasets processed through scAMF (Fig. 4B). This enhancement in neighborhood purity provides strong evidence of scAMF’s effectiveness in improving class distinction within the scRNA-seq data.

scAMF Provides Better Visualization.

In this part, we show that scAMF’s visualization performance surpasses those of widely used methods like T-SNE (28), uniform manifold approximation and projection (UMAP)(32), and the classical principle component analysis (PCA)(33). We present visualizations of two typical data, Kolodziejczyk and Usoskin, where different classes are indicated by various colors (Fig. 5A).

Fig. 5. 2D Visualizations and silhouette indexes of scAMF. (A) Color-coded representations of Kolodziejczyk and Usoskin data using T-SNE, UMAP, PCA, and scAMF. We also report the silhouette index (SI) for all plots, measuring the compactness of clusters within the same class and the separation between different classes. The T-SNE and UMAP results are presented under the optimal parameter settings corresponding to the highest SI. (B) A box plot of silhouette indexes for 8 clustering methods along with PCA, T-SNE, and UMAP across all datasets. scAMF significantly outperforms other methods in terms of silhouette index.

In analyzing the Kolodziejczyk data (Top row of Fig. 5A), PCA demonstrates limited class differentiation, with greater variation within classes than between them. T-SNE, akin to PCA, manages some separation but indicates potential within-class divisions, likely arising from batch effects in the dataset. UMAP, however, inaccurately splits the “2i” class into two distinct clusters. In contrast, scAMF effectively segregates the three classes, achieving clear distinction and creating densely packed clusters for each class, thus demonstrating its superior clustering capability.

The main challenge with the Usoskin data (Bottom row of Fig. 5A) is the overlapping nature of different classes, making complete separation challenging. Techniques like PCA, T-SNE, and UMAP struggle to distinguish between these classes adequately. T-SNE and UMAP show some improvement over PCA, yet this enhancement is still limited. In contrast, scAMF effectively differentiates all classes while minimizing the misplacement of cells. Notably, scAMF highlights the bridging role of the peptidergic nociceptors (PEP) group. scAMF’s visualization reveals a connection between PEP neurons and the non-peptidergic nociceptors (NP), tyrosine hydroxylase containing (TH), and neurofilament containing (NF) groups, while the latter three groups show distinct separation. This observation aligns with biological insights into PEP neurons (34). That is the marker genes in PEP neurons, including Tac1, Ntrk1, and Calca, may be highly expressed in other groups, while the expression of marker genes in the other three groups is relatively independent.

To assess the effectiveness of scAMF’s visualization capability, we have quantitatively compared it with other scRNA-seq analysis methods using the silhouette index. The silhouette index (35), ranging from −1 to 1, serves as a measure of visualization clarity and accuracy. Higher values indicate better visualization quality. This comparison is visually represented in a box plot (Fig. 5B), showing the silhouette index distribution across all datasets. The plot illustrates a notable shift of the scAMF box toward 1, indicating superior visualization performance compared to other methods. Detailed silhouette index scores for each method are available (SI Appendix, Table S5). Although scDHA may achieve a higher silhouette index in a few datasets, scAMF stands out for its consistent performance and lower variability, making it a more reliable choice for scRNA-seq data visualization. Further experiments have been conducted to verify how scAMF preserves global cell topology in comparison to PCA. Due to the transformation undergone by scAMF prior to manifold fitting, there has been a significant alteration in the neighborhood distribution. Consequently, the limited overlap with PCA is understandable, as depicted in SI Appendix, Fig. S3, where the median neighborhood overlap for most datasets does not exceed 50%. Furthermore, the exceptional performance of higher neighborhood class purity further underscores scAMF’s capability to discern nonlinear patterns effectively.

Supplementary experiments to optimize parameters for both T-SNE and UMAP are illustrated in SI Appendix, Fig. S4. Regardless of the parameter configurations, scAMF consistently demonstrates superior performance over T-SNE and UMAP in both SI and visualization outcomes.

Materials and Methods

Materials.

Our study collects 25 diverse scRNA datasets centered on cellular and molecular biology, encompassing embryonic, fetal, and adult tissues from humans and mice. These datasets offer valuable insights into various stages of development and disease, each with unique characteristics and sourced from different studies, as detailed by their respective Gene Expression Omnibus (GEO) accession codes and research citations. They provide a rich resource for understanding biological processes and heterogeneity at a single-cell level, employing a range of protocols, including Smart-Seq2, STRT-Seq, and 10X Genomics. All datasets are freely available for download, along with detailed information for each dataset (SI Appendix, F and Table S1).

Methods.

Suppose we have a dataset X={x(i)}i=1N, where X represents a collection of scRNA-seq raw data. Each vector x(i) corresponds to the expression values [x1(i),…,xD(i)] of the i-th cell across D genes. The transformed data are denoted as Y={y(i)}i=1N, where each transformed cell expression y(i) consists of [y1(i),…,yD(i)], still encompassing D genes. Finally, the fitted data are represented by Z={z(i)}i=1N.

Raw data transformation.

Raw data are transformed in parallel, allowing us to select the most suitable transformation. As we will show, it would achieve the best possible condition for downstream analysis. While multiple transformation methods can be chosen, we apply three of them in our experiments.

Value-to-rank transformation (T1). Instead of looking at the actual levels of gene expression, the value-to-rank transformation focuses on where each gene stands in the order within a cell. This method transforms the expression levels to their respective ranks within each cell expression vector. Specifically, the transformation is defined component-wisely with each yk(i) being the rank of xk(i) in the ascending order of x(i).

Unit-vector transformation (T2). To minimize the impact of varying expression magnitudes across cells, this transformation normalizes the expression vectors of each cell to have a unit Euclidean norm. For an expression vector x(i), the normalized vector is given by y(i)=x(i)/‖x(i)‖2 where ‖x(i)‖2 represents the Euclidean norm of x(i).

Logarithmic transformation (T3). Logarithmic transformation is also applied to rectify the skewness commonly present in scRNA-seq data. Each expression vector xk(i) in vector x(i) is transformed as follows yk(i)=log2(1+xk(i)).

This approach stabilizes variance across varying expression levels, making the data more suitable for subsequent analysis.

Manifold fitting.

The datasets transformed are labeled as Y1,⋯,YnT, respectively. For a transformed data Y∈Yi:i=1,...,nT, with D as the dimension of the ambient space (i.e., the number of genes), the manifold hypothesis assumes that the Y can be viewed as a noisy representation of the underlying low-dimensional manifold M of dimension d (d≪D). This implies that the true (noiseless) samples lie on a simpler, more basic structure but are obscured or altered due to noise. Hence, the fitting aims to construct an estimated manifold M^ from the noisy data Y, or equivalently, approximately project Y onto M^. Specifically, the main steps of the manifold fitting algorithm (23) include direction estimation and projection estimation.

Direction estimation.

The purpose of finding the projection direction is to determine the direction where the point should move closer to the underlying manifold. For any point y, its projection direction can be expressed as dy=F(y)−y, where F(y) is obtained through the local contraction of y. In ref. 23, F(y) is defined by the following formula:F(y)=∑iαi(y)y(i),

with αi(y) defined as:[1] α~i(y)=(1−‖y−y(i)‖22r02)k,‖y−y(i)‖2≤r0;0,otherwise;α~(y)=∑i∈Iyα~i(y),αi(y)=α~i(y)α~(y),

and k>2 (3 by default) ensures smoothness, and Iy contains indices of members in BD(y,r0). In our implementation, we refine F(y) by addressing two key aspects. 1) recognizing the limitations of Euclidean distance for scRNA-seq data, we have adopted shared nearest neighbor metrics—a method better suited for high-dimensional data—to accurately determine the neighboring points involved in computation F(y). The function of r0 is substituted by a parameter p, which is introduced in detail in the following session. 2) we have set the parameter to k=0 in Eq. 1 to simplify the calculating process. This simplification is practical, as our objective is solely to project the points without the need for the smoothness of the output manifold. We also explored the impact of smoothness on clustering outcomes, with findings presented in SI Appendix, Fig. S5. It’s observed that an increase in the value of k leads to diminished clustering performance. Specifically, for large k values (greater than 1), the weights of points ranked lower within B(i) exponentially decrease, minimally contributing to the fitting of F(y(i)). In other words, when k is large, the resulting smoothness means that not all neighbors selected by the refined p-nearest neighborhood contribute to the fitting of F(y(i)), and consequently impacts the clustering results. Hence, in scAMF, we opt to compromise on smoothness to ensure superior clustering performance. The rest is the detailed calculation of F(y).

– Computing the shared nearest neighborhood. For each y(i), we denote Np(i) as p (default by 15) nearest neighborhood of y(i), as determined by a given metric. For all datasets we examined, robustness to the parameter p in scAMF is tested. Detailed information is available in SI Appendix, Fig. S6. Np(i) contains the key sample points where y(i) should pay attention to. Hence, the shared nearest neighborhood (SNN) of y(i) and y(j) is defined as: SNNi,j=Np(i)∩Np(j).

– Determine the neighborhood. Using the SNN of y(i) and y(j), the refined p-nearest neighborhood of y(i) is defined by B(i)=argmaxS⊂Y,|S|=p∑y(j)∈SSNN(i,j).

– Calculate the projection direction. The projection direction F(y(i)) of y(i) is defined by F(y(i))=1|B(i)|∑y(j)∈B(i)y(j).

Projection estimation.

Upon determining F(y(i)), we construct a cylinder centered at y(i) with F(y(i))−y(i) as its major axis, following ref. 23. The weighted average of the samples within this cylinder, denoted as G(y(i)), serves as an estimate for the projection. To accommodate the ultrahigh nature and irregular noise patterns in scRNA-seq data, we simplify G(y(i)) by identifying the point of maximum density along the line connecting y(i) and F(y(i)) and using this point as a substitute for the original G(y(i)). Precisely, the simplification isG(y(i))=argmaxytρ(yt),

whereyt=y(i)+t(F(y(i))−y(i))

andρ(y(i))=1∑y(j)∈B(i)‖y(i)−y(j)‖22

quantifies the density of y(i), with higher values indicating denser regions. This modification simplifies the analysis while accommodating the complex nature of scRNA-seq. The fitted data are then represented byZ={z(i)}i=1N,z(i)=G(y(i)).

Unsupervised clustering and validation.

Upon outputting Z, one may simply utilize fairly basic clustering methods. These could include but are not limited to

Spectral Clustering (C1): group the data by analyzing the connectivity (adjacent matrix) among data points. Subsequently, it uses the k-means algorithm to perform clustering within this reduced-dimensional space.

Agglomerative Hierarchical Clustering (C2–C4): creates a hierarchy of clusters, with each being characterized by its approach to linkage. Single linkage calculates the distance between clusters based on the closest pair of points, one from each cluster. Average linkage takes the average distance between all pairs of points across two clusters. Centroid linkage defines the distance between clusters as the distance between their respective centroids.

We remark that the clustering module in scAMF necessitates a predetermined number of cell types, denoted as K. Researchers can provide this number based on biological insights. To ensure a fair comparison, we employ the true number of cell types for clustering in all compared methods that require K as input (SI Appendix, D).

By implementing various transformations Ti and clustering methods Cj, where 1≤i≤nT and 1≤j≤nC, multiple distinct clustering outcomes Cij can be obtained. Each outcome is derived from the fitted data Zi and categorizes cells into K groups. To report the optimal clustering result, we further introduce an UVI based on the similarity ratio to measure the justifiability of each Cij. We denote the label of the h-th cell in Cij by Cij(h). The similarity ratio estimates the clustering quality based on two primary factors: average intraclass and interclass similarity. Specifically, the similarity between z(h) and z(ℓ), denoted as s(h,ℓ), is defined as:s(h,ℓ)=1ifz(h)∈Nq(ℓ)orz(ℓ)∈Nq(h);0otherwise,

where q is adaptively determined as min{N/10,40}. The average intraclass similarity, represented as a(z(h),Cij), evaluates the cohesion within a cluster asa(z(h),Cij)=1|{(h,ℓ):Cij(h)=Cij(ℓ)}|∑Cij(h)=Cij(ℓ)s(h,ℓ).

Furthermore, the interclass similarity, measuring the separation between different clusters, is denoted byb(z(h),Cij)=1|{(h,ℓ):Cij(h)≠Cij(ℓ)}|∑Cij(h)≠Cij(ℓ)s(h,ℓ).

Hence, the UVI through the similarity ratio is defined byUVI(Z,Cij)=∑ha(z(h),Cij)∑hb(z(h),Cij).

A higher UVI indicates superior clustering performance, marked by greater similarity within clusters than between them. We select the clustering result with the highest UVI from all Cij, and denote it by C, as the final output (refer to Algorithm 1).

Discussion

In this study, we present scAMF, a framework enhancing scRNA-seq analysis through manifold fitting. scAMF insightfully represents the hidden low-dimensional structure, addressing technical variability and biological noise inherent in scRNA-seq data. Rigorous experimentation demonstrates the superior performance of scAMF in recovering distorted RNA expression data and improving clustering compared to existing techniques.

In our experiment, the true number of cell types is used as input to implement and evaluate scAMF and other clustering methods. Determining the number of cell types is an open problem in RNA-seq analysis. There are two strategies for determining the number of cell types. The first strategy utilizes statistical techniques like the silhouette coefficient (35) or Dunn index (36), which may sometimes overestimate the true number. Alternatively, as suggested in ref. 37, biologists should predetermine the number of cell types based on their research needs, opting for more clusters for detailed exploration or fewer for simplified analysis. While biologist-provided counts might slightly deviate from the true number of cell types, scAMF exhibits considerable flexibility, effectively managing both overestimation and underestimation (SI Appendix, Table S7). For instance, even when the number of cell types is overestimated by three, scAMF still achieves a notable average ARI of 0.66 for all datasets, outperforming other methods using the true number of cell types.

Future research directions for the manifold fitting approach used in scAMF could focus on expanding its application to other omics data like proteomics and metabolomics to tackle similar high-dimensionality and noise challenges. Additionally, enhancing the computational efficiency of scAMF, particularly for processing increasingly large datasets, could be achieved through more efficient algorithms or parallel computing. Moreover, integrating scAMF with emerging spatial transcriptomics data would offer a more comprehensive understanding of tissue architecture and cell-to-cell interactions, enabling spatial mapping of cell types and states.

Supplementary Material

Appendix 01 (PDF)

Z.Y. has been supported by Singapore Ministry of Education Tier 2 grant (A-0008520-00-00 and A-8001562-00-00) and Tier 1 grant (A-0004809-00-00 and A8000987-00-00) at the National University of Singapore; Y.L. is supported by Research Scholarship at the National University of Singapore. Z.Y. thanks and appreciates the constructive discussion with Jun Liu at Harvard University.

Author contributions

Z.Y. and S.-T.Y. designed research; Z.Y., B.L., and Y.L. performed research; Z.Y., B.L., and Y.L. contributed new reagents/analytic tools; Z.Y., B.L., and Y.L. analyzed data; and Z.Y., B.L., Y.L., and S.-T.Y. wrote the paper.

Competing interests

The authors declare no competing interest.

Data, Materials, and Software Availability

The MATLAB-based implementation of scAMF, encompassing a demo of the scAMF pipeline, the three datasets used in our study, their corresponding results, and all the evaluation functions employed, is readily available for download at https://github.com/zhigang-yao/scAMF (24). All other data are included in the manuscript and/or SI Appendix.

Supporting Information

Reviewers: S.H., Georg August University of Göttingen; and J.J.L., University of California, Los Angeles.
==== Refs
1 F. Tang , mRNA-Seq whole-transcriptome analysis of a single cell. Nat. Methods 6 , 377–382 (2009).19349980
2 R. Grindberg , RNA-sequencing from single nuclei. Proc. Natl. Acad. Sci. U.S.A. 110 , 19802–19807 (2013).24248345
3 Y. Hao , Integrated analysis of multimodal single-cell data. Cell 184 , 3573–3587 (2021).34062119
4 F. Wolf, P. Angerer, F. Theis, SCANPY: Large-scale single-cell gene expression data analysis. Genome Biol. 19 , 1–5 (2018).29301551
5 S. Sachs , Targeted pharmacological therapy restores β-cell function for diabetes remission. Nat. Metab. 2 , 192–209 (2020).32694693
6 H. Mathys , Single-cell transcriptomic analysis of Alzheimer’s disease. Nature 570 , 332–337 (2019).31042697
7 Y. Katzenelenbogen , Coupled scRNA-seq and intracellular protein activity reveal an immunosuppressive role of TREM2 in cancer. Cell 182 , 872–885 (2020).32783915
8 Z. Miao, B. Humphreys, A. McMahon, J. Kim, Multi-omics integration in the age of million single-cell data. Nat. Rev. Nephrol. 17 , 710–724 (2021).34417589
9 L. Moses, L. Pachter, Museum of spatial transcriptomics. Nat. Methods 19 , 534–546 (2022).35273392
10 A. Regev , The human cell atlas. eLife 6 , e27041 (2017).29206104
11 P. Lin, M. Troup, J. Ho, CIDR: Ultrafast and accurate clustering through imputation for single-cell RNA-seq data. Genome Biol. 18 , 1–11 (2017).28077169
12 D. van Dijk et al., MAGIC: A diffusion-based imputation method reveals gene-gene interactions in single-cell RNA-sequencing data. bioRxiv [Preprint] (2017). 10.1101/111591 (Accessed 7 March 2023).
13 M. Huang , SAVER: Gene expression recovery for single-cell RNA sequencing. Nat. Methods 15 , 539–542 (2018).29941873
14 W. Li, J. Li, An accurate and robust imputation method scImpute for single-cell RNA-seq data. Nat. Commun. 9 , 997 (2018).29520097
15 T. Stuart , Comprehensive integration of single-cell data. Cell 177 , 1888–1902 (2019).31178118
16 V. Kiselev , SC3: Consensus clustering of single-cell RNA-seq data. Nat. Methods 14 , 483–486 (2017).28346451
17 V. Blondel, J. Guillaume, R. Lambiotte, E. Lefebvre, Fast unfolding of communities in large networks. J. Stat. Mech. Theory Exp. 2008 , P10008 (2008).
18 D. Tran , Fast and precise single-cell data analysis using a hierarchical autoencoder. Nat. Commun. 12 , 1029 (2021).33589635
19 X. Li , Deep learning enables accurate clustering with batch effect removal in single-cell RNA-seq analysis. Nat. Commun. 11 , 2338 (2020).32393754
20 K. Mohammed, H. Narayanan, Manifold learning using kernel density estimation and local principal components analysis. arXiv [Preprint] (2017). https://arxiv.org/abs/1709.03615 (Accessed on 20 March 2023).
21 C. Fefferman, S. Ivanov, Y. Kurylev, M. Lassas, H. Narayanan, “Fitting a putative manifold to noisy data” in Conference on Learning Theory, B. Sébastien, P. Vianney, R. Philippe, Eds. (PMLR, 2018), pp. 688–720.
22 C. Fefferman, S. Ivanov, M. Lassas, H. Narayanan, Fitting a manifold of large reach to noisy data. arXiv [Preprint] (2021). https://arxiv.org/abs/1910.05084 (Accessed 20 March 2023).
23 Z. Yao, J. Su, B. Li, S.-T. Yau, Manifold fitting. arXiv [Preprint] (2023). https://arxiv.org/abs/2304.07680 (Accessed 22 March 2023).
24 Z. Yao, B. Li, Y. Lu, S.-T. Yau, zhigang-yao/scAMF. GitHub. https://github.com/zhigang-yao/scAMF. Deposited 21 May 2024.
25 Z. Yao, J. Su, S.-T. Yau, Manifold fitting with CycleGAN. Proc. Natl. Acad. Sci. U.S.A. 5 , e2311436121 (2023).
26 H. Wiechers, B. Eltzner, K. V. Mardia, S. F. Huckemann, Learning torus PCA-based classification for multiscale RNA correction with application to SARS-CoV-2. J. Roy. Stat. Soc. Ser. C Appl. Stat. 72 , 271–293 (2023).
27 A. Kolodziejczyk , Single cell RNA-sequencing of pluripotent states unlocks modular transcriptional variation. Cell Stem Cell 17 , 471–485 (2015).26431182
28 L. Van der Maaten, G. Hinton, Visualizing data using T-SNE. J. Mach. Learn. Res. 9 , 2579–2605 (2008).
29 L. Ertoz, M. Steinbach, V. Kumar, “A new shared nearest neighbor clustering algorithm and its applications” in Workshop on Clustering High Dimensional Data and its Applications at 2nd SIAM International Conference on Data Mining (2002), pp. 105–115.
30 M. Guo, S. Hui Wang, J. W. Potter, X. Yan, SINCERA: A pipeline for single-cell RNA-seq profiling analysis. PLoS Comput. Biol. 11 , e1004575 (2015).26600239
31 Y. Hao et al., Dictionary learning for integrative, multimodal and scalable single-cell analysis. Nat. Biotechnol. 41, 1270–1281 (2023).
32 L. McInnes, J. Healy, J. Melville, UMAP: Uniform manifold approximation and projection for dimension reduction. arXiv [Preprint] (2018). https://arxiv.org/abs/1802.03426 (Accessed 13 June 2023).
33 H. Hotelling, Analysis of a complex of statistical variables into principal components. J. Educ. Psychol. 24 , 417 (1933).
34 D. Usoskin , Unbiased classification of sensory neuron types by large-scale single-cell RNA sequencing. Nat. Neurosci. 18 , 145–153 (2015).25420068
35 P. Rousseeuw, Silhouettes: A graphical aid to the interpretation and validation of cluster analysis. J. Comput. Appl. Math. 20 , 53–65 (1987).
36 J. C. Dunn, A fuzzy relative of the ISODATA process and its use in detecting compact well-separated clusters. J. Cybern. 3 , 32–57 (1973).
37 L. Chen, W. Wang, Y. Zhai, M. Deng, Deep soft k-means clustering with self-training for single-cell RNA sequence data. NAR Genomics Bioinforma. 2 , lqaa039 (2020).
