Advertisement
Methods
Transparent Process
Open Access

A scalable approach to topic modelling in single-cell data by approximate pseudobulk projection

Sishir Subedi, Tomokazu S Sumida, Yongjin P Park  Correspondence email
Published 6 August 2024. DOI: 10.26508/lsa.202402713

Abstract

Probabilistic topic modelling has become essential in many types of single-cell data analysis. Based on probabilistic topic assignments in each cell, we identify the latent representation of cellular states. A dictionary matrix, consisting of topic-specific gene frequency vectors, provides interpretable bases to be compared with known cell type–specific marker genes and other pathway annotations. However, fitting a topic model on a large number of cells would require heavy computational resources–specialized computing units, computing time and memory. Here, we present a scalable approximation method customized for single-cell RNA-seq data analysis, termed ASAP, short for Annotating a Single-cell data matrix by Approximate Pseudobulk estimation. Our approach is more accurate than existing methods but requires orders of magnitude less computing time, leaving much lower memory consumption. We also show that our approach is widely applicable for atlas-scale data analysis; our method seamlessly integrates single-cell and bulk data in joint analysis, not requiring additional preprocessing or feature selection steps.

Introduction

Background

High-throughput single-cell sequencing has gained popularity and has been successfully applied to many recent studies to understand cellular contexts. Not only the number of studies that involve single-cell sequencing but also the number of cells that a single study can profile has dramatically increased from 100s to millions. Advancements in computational methods have necessarily accompanied such technological advances in order to provide necessary operations in almost every step of data analysis—from data procurement to final statistical inference (Kharchenko, 2021). However, typical computing infrastructures are often not well equipped with powerful computing units and high-capacity memory despite the fact that many advancements in computational methods assume such large resources are readily accessible (Heumos et al, 2023; Zhang et al, 2024). Training a sophisticated model on millions of single-cell data vectors, using a stochastic gradient descent algorithm, not surprisingly demands much more resources than fitting a model on data with 1,000s of cells and a 1,000 features.

After the quantification of molecular features within each cell, a conventional single-cell analysis (Stoeckius et al, 2017; Stuart & Satija, 2019; Hao et al, 2021) generally turns to dimensionality reduction, such as principal component analysis (PCA) (Hotelling, 1933; Jolliffe, 1986), nearest neighbour graphs were constructed by matching cells based on the reduced data, and cell clusters are often resolved by a graph-based clustering method, such as the Louvain (Blondel et al, 2008) or the Leiden algorithm (Waltman and van Eck, 2013; Traag et al, 2018 Preprint). The complexity of an exact PCA algorithm (singular value decomposition) scales linearly with the number of genes/features and quadratically with the number of cells (Tsuyuzaki et al, 2020), designing a scalable latent data modelling is crucially important (Sun et al, 2019). Should a large number of cells be analysed, finding hidden gene regulatory programs/factors/topics is a computationally intensive task. It is customary to conduct an unvetted feature selection for highly variable genes to reduce memory footprint. By doing so, we critically assume that biologically relevant signals are most pronounced in data and further take that the rank of gene expression variance of the selected genes is somewhat stable across different data sets.

Related work

How can we uncover latent gene programs without compromising the intrinsic dimensionality? Several approaches have been suggested to overcome the scalability issues of the embedding or dimensionality reduction step before or after building neighbour graphs, and the ultimate outcome is to resolve clusters of cells using a graph-based clustering method. Most existing methods address the question of how we select fewer representative cells to ease the computational burden of downstream analysis, not sacrificing the accuracy of the final clustering results (Hie et al, 2019b; Abdelaal et al, 2020; Dhapola et al, 2022; Wei et al, 2022). A geometric sketching method was first coined to randomly sample cells while uniformly covering an underlying and unknown metric space (Hie et al, 2019b). A similar idea was explored in natural language processing but more focused on defining the convex hull (boundary) of word occurrence topics over document corpus (Mimno & Lee, 2014). More recent scalable approaches developed by the genomics community also build on the same premise that we can approximately capture the overall clustering patterns across millions of cells by choosing representative (anchor) cells wisely (Abdelaal et al, 2020; Dhapola et al, 2022; Wei et al, 2022). After building nearest neighbour graphs of the anchor cells, the rest of the algorithmic steps in all these approaches largely focus on adding new cells (vertices) to the original nearest neighbour graphs of the initial anchor cells and performing label propagation of cell-type annotations.

Our approach

Instead, here, we focus on developing a scalable approach for probabilistic topic modelling. We deal with high dimensionality and large sample size by mapping high-dimensional expression/activity vectors onto the lower dimensional topic space via random projection. Probabilistic topic modelling decomposes a large count data matrix into an interpretable dictionary matrix (topic-specific activities across 10s of 1,000s of genes/features) and a topic proportion matrix (each cell’s attribution to the topics) (Blei et al, 2003; Dieng et al, 2020). Large-scale single-cell analysis can be made quite rapidly by considering the data generation process from a rather different perspective. Instead of designing a new type of complex model while leaving its statistical inference up to machine-learning libraries, we propose a new framework in which a fundamental topic model can be estimated in an accurate, robust, and scalable way.

Why topic modelling

Fitting a probabilistic topic model has long been considered a principled and intuitive approach to uncovering patterns hidden underneath count data derived from high-throughput sequencing. A grade of membership model (or admixture) was first coined in genetics while trying to identify population structures manifested in genetics variant counts across the genome (Pritchard et al, 2000; Novembre, 2016). In the same way, multi-tissue, multi-individual gene expression patterns were analysed by the same type of model (Dey et al, 2017). More recently, the Poisson matrix factorization (PMF) was shown to be an equivalent problem, and a more efficient method based on alternating regression estimation was suggested (Carbonetto et al, 2021 Preprint, 2023 Preprint). In single-cell genomics, especially for sparse DNA accessibility data, embedding methods based on topic modelling (or latent semantic indexing), such as ArchR (Granja et al, 2021) and cisTopic (Bravo González-Blas et al, 2019), were ranked in the top lists of a recent benchmark study (Chen et al, 2019). A topic model approach based on a deep variational autoencoder model, called an embedded topic model, or ETM (Dieng et al, 2020), has been successfully used in single-cell RNA-seq modelling (Zhao et al, 2021b; Subedi & Park, 2023; Zhang et al, 2023).

Our work

We present a new scalable and versatile method that can quickly recover cellular topics with modest computing resources. We term our method ASAP, short for Annotating a Single-cell data matrix by Approximate Pseudobulk estimation. Several existing works and observations inspired us. First, repeatedly applying random projection operation onto high-dimensional feature vectors may result in quick and moderately accurate cell clustering patterns (Wan et al, 2020). Second, pseudobulk data derived from aggregating within cell types and states often behave similar to typical RNA-seq bulk data. Third, a grade of membership, or an equivalent non-negative matrix factorization (NMF), method is powerful enough to dissect cellular topics from bulk sequencing profiles (Dey et al, 2017) and single-cell sequencing data (Carbonetto et al, 2021 Preprint).

Results

Overview of ASAP

Briefly, ASAP will handle a matrix factorization problem of massive single-cell data in the following three steps: (1) we randomly project cells onto the low-dimensional space (Fig 1A) to sort them through a binary classification tree (Fig 1B), depending on the sign of the random projection (RP) values, and construct a pseudobulk (PB) data matrix by aggregating cells landed in the same termini of the binary sorting tree. As suggested by the previous work (Wan et al, 2020) and demonstrated in the UMAP (Becht et al, 2018) (Uniform Manifold Approximation and Projection), simple RP operations can help cells with a similar expression variation group together (Fig 1B). It is also encouraging to find that the resulting PB data already show that cells of similar cell types are naturally enriched within the same PB sample (the structure plot of Fig 1B, cell-type proportion by PB samples). (2) Followed by the NMF, or, more precisely, the PMF, we can decompose the PB expression data (Ypb) into topic-specific gene expression (dictionary) matrix (β) and topic loadings/proportions for each PB sample (θpb). The non-negativity constraint is natural to the gene expression count data, and the additivity of factors generates biologically interpretable solutions. (3) Treating the dictionary matrix as a design matrix in a non-linear, non-negative regression problem, we can quickly deconvolve cell-level topic proportions (Fig 1D). The resulting topic proportion matrix often results in markedly improved UMAP results (Fig 1B versus Fig 1D). Moreover, the dictionary matrix β and two topic proportion matrices θ, estimated by the PMF algorithm (Gopalan et al, 2014, 2016; Levitin et al, 2019), clearly exhibit modular structures in both sides (see the Materials and Methods section for technical details).

Figure 1.
Figure 1. Schematic overview of the ASAP framework.

(A, C, D) ASAP proceeds in three steps: collapsing cell count data into a manageable pseudobulk matrix (A), factorizing the pseudobulk matrix by the Poisson matrix factorization (C), and regressing cell-level data onto the topic space (D). For illustration purpose, we show intermediate results in processing pancreatic ductal adenocarcinoma data, consisting of 9,887 cells (Zhao et al, 2021a). (A) Raw gene expression matrix can be projected on the r-dimensional space of RP variables (right) by multiplying the observed cell count matrix (middle) with some univariate Gaussian matrix (left). (B) According to the sign of each random projection (RP) direction, switching to the left and right at each level, cells are sorted via a perfect binary tree with the depth =r (RP variables) and aggregated into pseudobulk samples. (C) Non-negative matrix factorization method decomposes the pseudobulk data matrix into a dictionary (gene by topic/factor) and factor loading matrix (sample by topic/factor). We represent each pseudobulk sample as a non-negative linear combination of non-negative topic-specific gene expression topics. (D) Treating the dictionary matrix as a design matrix, ASAP revisits the original gene expression vectors of many cells and regresses them on the topic space with the same non-negative constraints on the regression coefficients, which can be used to identify cell-level clustering patterns (Uniform Manifold Approximation and Projection on the right).

ASAP accurately estimates cell clusters

In single-cell data analysis, it is implicitly assumed that we collect a sufficient number of cells within each cell type, and these cells show more similar expression patterns within the same cell type than between different cell types. Some intrinsic, biological similarity measures can group cells of identical cell types/states. Had dimensionality reduction or latent embedding methods well preserved cells’ biological similarity, clustering results based on the latent factors/topics should have been closely matched with the groups based on ground truth cell types. Hence, we evaluate the performance of different embedding methods in terms of the correspondence between cell groups found by clustering each method’s latent factors/topics and ground truth cell types.

We generated realistic simulation data treating cell type–specific gene expression profiles in cell sorting experiments (Schmiedel et al, 2018). We used the sorted bulk sequencing data instead of other single-cell data sets that are computationally annotated by some clustering methods. To gain insights into variability induced by sequencing depths and sample size, we varied the number of cells from 1,000 to 10,000, while also varying the proportion of cell type–specific gene expression and background variation (see the Materials and Methods section). Because we are interested in comparing the estimated set membership (based on clustering) with actual cell-type membership, we used three performance metrics widely used to evaluate clustering tasks: adjusted Rand index (Rand, 1971), normalized mutual information (Vinh et al, 2010), and purity score (Liu et al, 2020; Kriebel & Welch, 2022).

Random projection well preserves underlying cell-type identities

We first wanted to understand the general behaviour of random projection operation. A similar analysis was conducted in previous studies, such as Wan et al (2020). As we increase the purity (ρ) level (Fig 2A), patterns emerge, and cells form groups (Fig 2B). Because visual inspection via UMAP is not rigorous enough, we evaluated the accuracy of the Leiden clustering results (Traag et al, 2018 Preprint) based on different configurations of latent factors derived from PCA and RP with a different number of components/factors (Fig 2E). Although the performance of RP data generally lags behind that of PCA by significant statistical margins, noting that PCA is more computationally intensive, we found that RP operation can be used as a good-enough-approximate method for the subsequent steps. Conversely, it also shows that RP alone is insufficient to capture a detailed view of cellular diversity, but the results of RP need to be refined.

Figure 2.
Figure 2. ASAP accurately estimates cell clusters.

Extensive benchmark experiments confirm that our approximating method performs better or at least as well as existing methods. (A) We generated realistic simulation data treating cell type–specific gene expression profiles in cell sorting experiments (Schmiedel et al, 2018) as a gold standard while applying a different noise level, or conversely increasing the proportion of cell type–specific variance from ρ = .01 to ρ = 1 (the top: ρ = 0.01; the middle: ρ = 0.5; the bottom: ρ = 1). (B) We can reduce dimension by multiplying the simulation count matrix with a random Gaussian matrix (see Fig 1A and Materials and Methods section). Uniform Manifold Approximation and Projection visualization shows how cell type–specific clusters emerge as we increase the proportion of cell-type variance from the top (ρ = 0.01) to bottom (ρ = 1). Each dot represents a cell with its originating cell type (colour). (C) Poisson matrix factorization decomposes a pseudobulk data matrix derived from the previous random projection step, yielding a topic-specific gene/feature frequency matrix (β) and topic proportions of each sample (θ). Here, we show the topic proportion matrices (θ), where each colour bar (x-axis) represents each sample’s probabilities of topics (y-axis; colours), and they sum to one. (D) Cell-level topic regression refines cell type–specific clusters. We show the Uniform Manifold Approximation and Projection based on the final topic regression results. (E) We tested the impact of different choices of RP dimensions on cell-type clustering performance. We compared clustering results derived from different embedding methods with the gold standard cell-type annotations (cell-sorted data). We report three different metrics on the y-axes—ARI (adjusted Rand index, the left), NMI (normalized mutual information, the middle), and purity (the right) while varying the noise levels. For all the metrics, the higher, the better. RP-k: the Leiden clustering of cell–cell 15-nearest neighbour graphs built on k-dimensional random projection results. PC-k: the Leiden clustering of cell–cell 15-nearest neighbour graphs built on top k principal components. (F) We compared the accuracy of ASAP’s final topic embedding results against other non-negative matrix factorization (NMF)–based methods, such as ASAP-full (the Poisson factorization directly applied to full data), Liger (Welch et al, 2019; Kriebel & Welch, 2022), NMF function, sklearn.decomposition.NMF, implemented in the scikit-learn library (Pedregosa et al, 2011), and a standard clustering pipeline implemented in the Scanpy library (Wolf et al, 2018). X-axis: the proportion of cell-type variance; y-axes: ARI (adjusted Rand index, the left), NMI (normalized mutual information, the middle), and purity (the right); for all the metrics, the higher, the better.

PMF methods achieve a higher level of robustness

Pseudobulk data

As a refinement step, we used the PMF (Gopalan et al, 2014, 2016; Levitin et al, 2019) applied to the PB data matrix (ASAP) or the full data matrix (ASAP-full/PMF). The PMF will decompose the resulting non-negative PB matrix into topic signatures (β) and topic proportion profiles (θ, Fig 2C). So, it was natural to compare with other NMF-based methods (Fig 2F), which include Liger (Welch et al, 2019; Kriebel & Welch, 2022), a standard NMF method (sklearn.decomposition.NMF) in the scikit-learn library (Pedregosa et al, 2011), along with a conventional Leiden clustering pipeline implemented in the Scanpy library (Wolf et al, 2018). All the methods were summarized into the clustering results, and cell memberships were tested against the gold standard (see the Materials and Methods section for details). Conceptually, we expect good clustering results distinctive clouds of points/cells in the topic space as visualized in UMAP (Fig 2D). The ASAP, Liger, and Scanpy’s clustering method naturally incorporates the Leiden clustering as a final step. For the Leiden-based methods, we varied the resolution parameters (0.1, 0.25, 0.5, 0.75, and 1.0); our results were, albeit, invariant to the fine-tuning of the Leiden algorithm. For the scikit-learn NMF, we simply used the standard k-means algorithm implemented in the same library.

We found that the PB-based ASAP and ASAP-full clearly outperformed other NMF and clustering methods with statistically significant margins, especially in settings with moderate and high noise levels (Fig 2F), suggesting that the PMF algorithm, which both ASAP and ASAP-full were built on, is more robust than other NMF algorithms. The PMF is a Bayesian hierarchical model, resolves its parameters by posterior inference, and often yields more robust results than optimization-based algorithms. Perhaps a more important aspect of this benchmark result is that our ASAP method based on PB data is as competent as the PMF fitted on full data (the red and orange curves in Fig 2F).

ASAP is a scalable method with little memory footprint

Along with the performance metrics (Fig 2), we also measured the runtime (in minutes) and memory footprint (in gigabytes/GB) of each method in the same workstation. We allocated the same resources for all the experiments, fixing the available memory to 32 GB and eight CPU cores. As long as the allocated memory permits, we varied the number of cells from 1,000 to 100,000. In the first batch of experiments, varying the number of cells from 1,000 to 10,000, ASAP-full (vanilla PMF) often took the most time (Fig 3A) and emended high memory to keep track of multiple copies of variational parameters (Fig 3B), followed by Liger and scikit-learn NMF methods. However, ASAP was consistently fast enough to get the results within a minute and kept the peak memory usage under 2 GB until the 10k cell experiments.

Figure 3.
Figure 3. ASAP achieves scalability in terms of computation time and memory footprint.

(A, B) Plot shows (A) the elapsed time (in minutes) and (B) peak memory usage (in GB) by non-negative matrix factorization–based methods with different numbers of cells ranging from 1,000 to 100,000 as a function of the number of cells. For the cases with the number of cells exceeding 10,000, we focused on comparisons between ASAP and Liger because our computing resource was insufficient to conduct full matrix factorization with a full set of features using non-negative matrix factorization implemented in the scikit-learn (Pedregosa et al, 2011) and the Poisson matrix factorization (Gopalan et al, 2014, 2016; Levitin et al, 2019) (ASAP-full).

In the second batch, we focused on Liger and ASAP models, varying the number of cells from 20,000 to 100,000. Here, ASAP showed significantly faster runtime and lower memory usage. Liger already peaked its memory usage near 32 GB with 100,000 cells. Therefore, it was not included in further experiments. We also conducted further experiments with ASAP alone to measure runtime and memory usage and found them to scale linearly with the number of cells (Supplementary Figure). Regretfully, our computing resources were limited; we had to focus on comparisons between the ASAP and Liger for experiments with the number of cells exceeding 10,000. A full matrix factorization with a full set of features using NMF implemented in the scikit-learn (Pedregosa et al, 2011) and PMF (ASAP-full) required substantially large memory (>128 GB) while taking days of computing time.

ASAP recapitulates known cell type–specific gene programs in real-world data sets

We applied ASAP to three different large-scale scRNA-seq data sets to demonstrate that our method is scalable and accurate enough to recapitulate well-established cell type–specific gene expression patterns. We considered the following data sets as representative examples: (1) PBMCs (Hao et al, 2021) profiled on 161,764 cells and 20,729 genes (Fig 4A–C); (2) breast cancer mammary tissues (Wu et al, 2021) of 20,265 genes expressed on 100,064 cells (Fig 4D and E); and (3) the first draft of the Tabula Sapiens project (Tabula Sapiens Consortium et al, 2022) (Fig 4F and G), consisting of 483,152 cells and 58,604 genes/transcripts.

Figure 4.
Figure 4. ASAP accurately annotates cell types in real-world large-scale data.

(A, B, C, D, E, F, G) Data sets: PBMC data (Hao et al, 2021) (A, B, C), breast cancer mammary tissues (Wu et al, 2021) (D, E), and Tabula Sapiens data (F, G) (Tabula Sapiens Consortium et al, 2022). (A, B, C) We show the three steps of ASAP applied to the PBMC data (A, B, C). (A) Step 1: Uniform Manifold Approximation and Projection (UMAP) drawn directly from the RP results. The colours represent different cell types previously annotated (level 1) by the original authors (Hao et al, 2021). (B) Step 2: pseudobulk data (top) along with the Poisson matrix factorization result—a dictionary matrix β (bottom). We compared the gene activities in the β matrix with the known cell-type marker genes (Franzén et al, 2019). (C) Step 3: fast cell-level regression analysis to recover topic proportions of all the cells. UMAP (top) drawn from the topic proportion results θ (bottom). (D) UMAP drawn directly from the RP results of the breast cancer data (Wu et al, 2021). (E) UMAP based on the final topic proportion results of ASAP. The colours represent different cell types based on the original study. Because of the space limit, we show the topic proportion matrix θ in Fig S1. (F) UMAP directly drawn from the RP results of the Tabula Sapiens data (Tabula Sapiens Consortium et al, 2022). (G) UMAP based on the final topic proportion results of ASAP. The colours represent different cell types based on the original study (the broad level). Because of the space limit, we show the topic proportion matrix θ in Fig S2.

ASAP’s dictionary matrix β (topic-specific gene expressions) strongly enriches the known cell-type marker genes, PanglaoDB (Franzén et al, 2019), within the same topics (bottom, Fig 4B); hence, the resulting cell topic proportions clearly differentiate cells based on the actual cell-type annotations (Fig 4C). We performed gene set enrichment analysis (GSEA) within each topic by fgsea (Korotkevich et al, 2021 Preprint), which provides a faster and more accurate approximation of the traditional GSEA method (Subramanian et al, 2005). Interestingly, the same level of resolution was not observed in the first RP step (Fig 4A). Many PB samples were rather heterogeneous, mixing multiple cell types (top, Fig 4B).

Considering that no label information was given to the factorization step, our intermediate results demonstrate that the non-negative matrix factorization step is indeed important, and the previous PB construction captures necessary information while compressing data complexity drastically. In the other real-world data analysis, we also observed that PB construction provided rough ideas about cell type–specific patterns (Fig 4D and F); the subsequent PMF steps refine topic-specific gene activities to result in high-resolution cell-type annotations in the end (Fig 4E and G). We also noted that the topics do not always agree with the cell types, which are often best characterized by cell type–specific surface markers and transcription factors. For instance, the topics of the breast cancer data may implicate finer resolution of cancer subtypes (Fig 4E), or the different cell types could be strongly affected by shared topic-specific genes (Fig S1). For a large study, such as Tabula Sapiens, cell-type annotations are often manually curated by different researchers with different domain knowledge, for example, dividing the jobs tissues-by-tissue. Nonetheless, we found the original annotations generally correspond to topics (Fig S2) and are well separated in the topic space (Fig 4G).

Figure S1.
Figure S1. Topic proportion matrix θ of the breast cancer mammary tissues (Wu et al, 2021).
Figure S2.
Figure S2. Topic proportion matrix θ of the Tabula Sapiens (Tabula Sapiens Consortium et al, 2022).

Moreover, we acknowledge that UMAP is perhaps not an ideal way to show a clear separation of distinct cell types manifested in the high-dimensional gene/feature space. Richer information is captured within the dictionary (β) and topic proportion (θ) matrices. We report supplementary figures to show the resulting topic proportion matrices (θ, cell by topic in Figs S1 and S2) and the dictionary matrices as the supplementary material.

Pseudobulk samples correspond to bulk expression samples

We had expected that PB samples would behave similar to cell type–sorted bulk samples. However, the previous real-world PB data showed quite a level of heterogeneity, mixing multiple cell types—especially those closely related cell types, such as CD4+ and CD8+ T cells. Such results also show that RP alone is insufficient to investigate fine-resolution cellular diversities.

To better understand the characteristics of PB samples, we investigated single-nucleus RNA-seq (snRNA-seq) (Eraslan et al, 2022) and bulk RNA-seq data (GTEx Consortium, 2020) profiled by Genotype-Tissue Expression (GTEx) consortium. This snRNA-seq data set provides gene expression vectors on 209,126 cells sampled from 16 individuals across 10 tissue types, for which the GTEx consortium also provides publicly accessible data (GTEx Consortium, 2020).

We first analysed the snRNA-seq matrix with ASAP and obtained its dictionary matrix β (gene by factor). As can be seen in the GSEA against PanglaoDB (Franzén et al, 2019) marker gene sets by the fgsea method (Subramanian et al, 2005; Korotkevich et al, 2021 Preprint), these topics are significantly associated with known cell types (bottom, Fig 5A). We also found that our unsupervised learning can recover the broad cell-type annotations called by the original study (Fig 5B). Interestingly, most cell types are shared across multiple tissues, yet there still exist several tissue-specific cell types, such as multiple epithelial cell types and myocytes, perhaps more specialized in tissue-specific environments.

Figure 5.
Figure 5. ASAP deconvolves cell type–specific gene programs in tissue-level data in joint analysis with snRNA-seq data.

(D, E) Data sets: Genotype-Tissue Expression project’s single-nucleus RNA-seq data (Eraslan et al, 2022) and GTEx project’s bulk RNA-seq data (GTEx Consortium, 2020) (D, E). (A) Topic-specific gene activity matrix β trained on the pseudobulk data constructed in the GTEx snRNA-seq data (Eraslan et al, 2022) shows strong enrichment of known cell types in PanglaoDB (Franzén et al, 2019). β: topic (rows) by gene (columns). (B) Patterns of ASAP’s topic proportion estimates θ agree with the original annotations. θ: cell (rows) by topic (columns). (C) Uniform Manifold Approximation and Projection embedding of GTEx bulk data, only including tissue types present in the snRNA-seq data. Before the Uniform Manifold Approximation and Projection, we projected the bulk data onto the snRNA-seq topic space. The colours represent different tissues of origin. (D) snRNA-seq pseudobulk samples (magenta) are overlaid on the same snRNA-seq topic space where the bulk samples (grey) are projected previously. (E) Cell-type proportions estimated by the joint analysis of GTEx snRNA-seq and bulk data. X-axis: cell-type fractions; y-axis: different tissue types; the numbers in the segments: fraction of cell types in percentage.

Using the dictionary matrix β, we projected the bulk RNA-seq data onto the snRNA-seq topic space (Fig 5C and D); then, we overlaid the snRNA-seq PB data on top of the same coordinates. GTEx tissues are generally well separated from one another (Fig 5C), suggesting different cell-type factions as a primary axis of variation. Moreover, we found that the PB samples overlap with these bulk samples to a large degree. Such good coverage compelled the follow-up experiment to assess whether the projected coordinates could be a basis for cell-type deconvolution step. Simply, we were able to query 100 neighbouring snRNA-seq cells for each bulk sample and count the frequency of each cell type within the neighbours.

The neighbouring cell-type fractions can be interpreted as approximate deconvolution profiles (Fig 5E). Unlike marker gene–based deconvolution methods, no prior knowledge was elicited, and reference cell-type profiles were adaptively selected for each bulk sample. Our results suggest several expected patterns. Substantially, higher fractions of myocytes were found in the sun-exposed side of skin and skeletal muscles (brown colours). A special type of alveolar epithelial cells was specifically found in the lung tissues. Immune cells are present in all the tissues but much less than other cell types. The two types of oesophagus tissues, histologically distinctive, are markedly different in terms of the cell-type fractions. We found a more diverse mixture in the oesophagus mucosa, containing a much higher fraction of fibroblasts and pericytes. We further investigated whether the cell-type proportion estimates induce visually separable embedding results on the UMAP space (Fig S3). It is clearly shown that major cell types are clearly distinctive from one another, forming disconnected clusters of cells of each type (the top panel). Subtypes were also visually distinctive but not as clear as the major categories (the bottom panel), implicating room for further improvement in future work.

Figure S3.
Figure S3. Cell-level Uniform Manifold Approximation and Projection visualizations of Genotype-Tissue Expression project’s single-nucleus RNA-seq data (Eraslan et al, 2022).

Sensitivity to the choice of hyper-parameters

Finally, we tested how the first step of our ASAP framework can be affected by the number of random projections, as the performance of the subsequent steps depends on it. We simulated single-cell data based on pbmc3k data available in the 10x Genomics website (see the Materials and Methods section for the detailed simulation scheme). Here, we fixed the sequencing depth to 10,000 and very little noise level (ρ = 0.025) and assumed the existing cell-type annotations as a gold standard. As expected, an insufficient number of random projection dimensions led to degraded performance (Fig 6A). However, we also found that the performance quickly recovered with sufficiently deep random projection trees above six, which leads to 26 = 64 pseudobulk samples, and no significant performance gain was observed beyond the optimal level. In practice, we found ASAP with a depth of 10 or above generally performs consistently well enough to construct an informative pseudobulk data matrix (≤1,024 samples).

Figure 6.
Figure 6. ASAP with a sufficiently deep binary–sorting tree performs robustly despite shallow sequencing depth.

(A) Varying the number of random projections, we evaluate the accuracy of the clustering patterns resulted from ASAP. X-axis: the depth of random binary trees; y-axis: the accuracy (NMI, ARI, purity). Each experiment was repeated five times. The vertical error bars represent one SD from the mean estimates. NMI, normalized mutual information (Vinh et al, 2010); ARI, adjusted Rand index (Rand, 1971); purity defined in references (Liu et al, 2020; Kriebel & Welch, 2022). For all the metrics, the higher, the better. (B) Clustering performance-varying sequencing depth parameters from 400 to 10,000 per cell. We simulated the data using 3,000 PBMCs while sampling five cells with the noise level ρ = 0.025, and resolved clustering by the Leiden method with the resolution parameter 0.1. For ASAP, we set the number of random projections to 10. X-axis: sequencing depth of the simulated data; y-axis: the accuracy (NMI, ARI, purity).

We also noted that the quality of single-cell data summarized by the coverage of non-zero features within each cell could be an important factor. Based on the same simulation scheme, we conducted rather small-scale experiments but with different levels of sequencing depths, that is, the number of reads mapped within each cell, varying from 400 to 10,000 (Fig 6B). To our surprise, the ASAP (PMF with pseudobulk) and PMF (ASAP-full) substantially outperformed across all the conditions. All the methods found a hard time picking up the right cell type–specific patterns with lots of missing values; hence, the scores with a shallow sequencing depth are significantly lower than those with a deep sequencing depth. The drop was rather steep in other NMF-related methods (Liger and scikit) and Scanpy’s Leiden clustering, whereas ASAP and PMF methods tend to maintain high accuracy even with little information.

Discussion

We propose a novel approximation method that can quickly identify topic structures embedded in large-scale single-cell data. Our method builds on the rationale that the distance between two high-dimensional data points can be preserved in the low-dimensional space constructed by random with a high probability (Johnson & Lindenstrauss, 1984; Bingham & Mannila, 2001; Blum, 2006; Dasgupta & Freund, 2008). Based on this foundation, a similar strategy has been used to speed up a clustering procedure (Wan et al, 2020); however, the quality of RP-based clustering analysis demands multiple random projections and clustering steps, and the resolution of clustering results is limited and unable to delineate a subtle difference between similar cell types; hence, the benefit of reduced computation time and peak memory can be diminished, and clustering results are often left less interpretable and require further post hoc analysis.

Although our approach begins with random projection steps, which had been previously explored for clustering, our key idea focuses on building a randomly sorted pseudobulk data matrix. We may consider randomly sorted pseudobulk data as a type of “meta”-cell. The potential of meta-cells has been explored in other studies, such as the SEACells algorithm (Persad et al, 2022 Preprint), which suggests efficient and accurate archetypal analysis. Unlike noisy (incomplete) cell–cell interaction patterns, meta-cells better represent intrinsic cellular states and complete random drop-out observations to help ascertain relationships between cis-regulatory regions and target genes. Our work is built on the same premise. Knowing that many cells are stochastic realizations of a “meta-cell,” as long as we can construct meaningful aggregates of cells, it is perhaps more sensible to run computational methods on the resulting meta-cell (pseudobulk) data.

In realistic benchmark analysis based on sorted bulk RNA-seq profiles, we demonstrate that learning a topic model on randomly sorted pseudobulk data outperforms unsupervised cell-type annotation tasks. Unlike clustering results, our topic model parameters are immediately interpretable and directly comparable against known marker gene databases. Furthermore, we show that our approach resembles a data-generating scheme of bulk RNA-seq data, and the inference results are easily transferable for deconvolution analysis. Especially for a data set with shallow sequencing depth, our method tends to show robust performance, which is highly desirable for large-scale single-cell data that contain a large fraction of missing observations.

Our approach can be extended in many different ways. We can intuitively apply RP operations to reduce the feature space for epigenomics, genomics, and other types of sparse matrix data. Although aggregating count data within a PB sample generally works robustly, we can combine batch effect adjustment and advanced depth normalization methods to further improve the quality of the resulting PB data. If needed, multi-omics data integration can take place at a PB level followed by nearest neighbour matching (Haghverdi et al, 2018; Hie et al, 2019a; Polański et al, 2020) or even optimal transport (Demetci et al, 2022). Moreover, in the case of heterogeneous multi-modal data integration tasks (Hristov et al, 2022), we expect matching results within the randomly projected space to strike a balance between runtime and theoretical guarantee.

Joint analysis in single-cell omics has gained more attraction; a recent method, such as multi-modal topic modelling (Kazwini & Sanguinetti, 2024), can benefit from our ASAP framework, speeding up Bayesian inference and achieving robustness in the presence of unwanted missing values. Any machine-learning algorithm can be used to learn a topic-specific feature frequency matrix. A deep learning–based generative model can be more efficiently trained by alternatively applying RP and PB steps; proper modelling that incorporates prior knowledge of cell type-specific marker genes and dependency structures will also substantially improve inference results (Elyanow et al, 2020; Townes & Engelhardt, 2023).

Materials and Methods

Probabilistic topic modelling by ASAP

Goal: topic modelling for single-cell data analysis

Single-cell transcriptomic profiling results in a count matrix X∈R≥0D×n of n cells where each column vector lies in the D-dimensional space. The dimensionality D corresponds to the number of features, such as genes and non-coding RNA molecules. Each element Xi (j ∈ [n]) quantifies how many reads were mapped onto a gene i within a cell j; hence, each value is non-negative by nature and generated by a counting process. Throughout this work, we will indicate a gene by i ∈ [D], denote a cell/sample by j, and use k ∈ [K] for a topic index. In total, we have D features and n cells, and the goal is to estimate K topics, (β1, … ,βK), and the corresponding topic proportion vector, (θ1, … ,θn), for each cell.

Given that each cell j, the column vector xj, sits on a topic space, k ∈ [K], we can assign a topic proportion of a topic k to the cell j, namely, θjk values, to quantify which and how much each topic affected each cell. Defining a matrix β with each element βik for an expected expression value for a gene g on a topic k, we assume that an observed gene expression matrix X was generated by a linear combination of the two matrices—β and θ. More precisely, we haveEXij|β,θ=∑k=1Kβikθjk,where the equality holds proportionally up to some scaling factor. The goal is to estimate the unknown gene topic (program) dictionary (β) and the topic proportion matrix (θ) provided with an observed gene expression matrix X and the fact that distinctive gene expression programs/topics exist and are expressed differently according to each cell’s topic representation.

PMF

The PMF model (Gopalan et al, 2014, 2016; Levitin et al, 2019) treats the same non-negative factorization problem as the Bayesian inference problem, assuming that each element of the data matrix, Xij, was generated by the following generative model:

1. Sample topic k-specific gene/feature expression: βik ∼ Gamma(a0, b0).

2. Sample each cell/sample’s activity/loading on a topic k: θjk ∼ Gamma(a0, b0).

3. Constitute the sampling rate λij for gene/feature i in each cell/sample j as a linear combination of the two parameters, namely, λij = Σk βik θjk.

4. Sample gene expression for a gene i in a sample j from the Poisson distribution: Xij ∼ Poisson(λij).

More precisely, we define the probability density function (PDF) of Gamma(a0, b0) aspλ|a,b=baΓaλa-1exp-bλ,and the PDF of Poisson(λ) aspX|λ=λXX!exp-λ.

Remark: The scale of the β and θ parameters is arbitrary, and they can adapt to each other, for example, βkθk⊤=(sk−1βk). The lack of identifiability in scaling factors brings about an interesting equivalence relationship between the PMF and the conventional multinomial topic model likelihood (Carbonetto et al, 2021 Preprint). It is straightforward to translate the PMF parameters into the conventional topic model parameters by adjusting the scales of β and θ parameters adaptively, such that Σi βik =1 and Σk θjk = 1.

Detailed steps of the ASAP algorithm

Step 1. Random projection (RP) and pseudobulk (PB) construction

A full D×n data matrix X can be sparse and too big to be read and densified in local memory. In the first step, each column/cell vector xj is compressed to a smaller, manageable d × n matrix Q by randomly projecting the total data matrix d times onto the d-dimensional space (d ≪ D). For each random projection (k = 1, … ,d), we use isotropic D-dimensional Gaussian vector rk and take the dot product with each cell vector, namely,qk⊤←[rk⊤x1,…,rk⊤xn].

To ensure that each RP vector is independent of the others, we perform (economical) singular value decomposition (SVD) of the Q matrix and construct orthogonal d × n RP matrix Q∼. After applying row-wise standardization, for each cell j ∈ [n], we can quickly convert q∼j≡ to a vector of d binary codes by setting each k-th coordinate to 1 if the Q∼kj is positive; otherwise, it is 0 (Fig 1B). This binary code vector can be used as a binary sorting tree, literally converting a binary number to a decimal one, with which we can quickly scatter cells into 2d buckets (leaf nodes). In each bucket/leaf node, we can incrementally aggregate expression vectors to a PB sample, and each PB vector yl is later normalized to have uniform sequencing depths (l = 1, … , 2d). For simplicity, we ignore empty buckets; hence, the resulting PB matrix can have at most 2d columns/samples.

• Construct a raw RP matrix Qd×n ← Rd×DXD×n.

• Orthogonalize the Q matrix and build Q∼ by SVD and standardization.

• Make a binary B matrix by setting each element.Bkj={1Q∼kj>00otherwise,where k = 1, … ,d and j = 1, … , n.

• Convert this B matrix to a leaf node (bottom bucket) membership matrix:Llj={1∑k=1d2k−1Bkj=l0otherwise,where k = 1, … , d and j = 1, … , n.

• Aggregate cells within each leaf node into a PB sample:yl=∑j:Llj=1xj

• Normalize yl ← yl/∥y∥1×104.

Step 2. Variational inference of the PMF model

Next, we decompose the d × L PB matrix Y by the PMF model (L ≤ 2d), modelling each PB data point as the Poisson distribution with the rate parameter that can be decomposed into dictionary and loading matrices (see the previous subsection for details). However, the log-likelihood,F=∑ijYijlog(∑k=1Kβikθjk)−∑ij∑k=1Kβikθjk+constant,is non-conjugate with the Gamma distributions of the β and θ parameters. Following the previous variational inference algorithm (Gopalan et al, 2014, 2016; Levitin et al, 2019), we need to introduce auxiliary variables to construct evidence lower bound (ELBO) under the mean-field variational approximation (using Jansen’s inequality):F≥∑ijkzijkYijlogβikθjkzijk−∑ijβikθjk+constant=ΔL.

This ELBO now permits the following closed-form update equations for the variational distributions.

• Inference of the local/sample-specific topic proportions:θjk|z,β∼Gammaθjk|a0+∑i=1Dyijzijk,b0+∑i=1DEqβik.

• Inference of the global/dictionary parameters:βik|z,θ∼Gammaβik|a0+∑j=1Lyijzijk,b0+∑j=1LEqθjk.

• Inference of the auxiliary variable:log zijk←Eq[logθjkβik∑kθjkβik]+constant,such that∑k=1Kzijk=1.

Step 3. Prediction using latent variable factorization

Finally, given that we have the expected dictionary matrix, namely, E[β] and E[logβ], we revisit full single-cell data and recover cell-level topic proportion vectors, namely, θj = for each cell j 1, … , n. We could use the same variational inference algorithm, only skipping the update for the β parameters, but it would involve full-scale inference of the auxiliary variables and not be scalable in practice. Instead, we form a different approximation of the ELBO and calibrate the θ parameters by solving a massive array of regression problems:L≥∑j,k∑iYijtjklogβikθjktjk−∑i,j,kβikθjk,where ∑ktjk=1 for all j. Letting ρjk=ΔE[tjk] under the same variational approximation, we can derive two closed-form update equations. Firstly, we havelog ρjk=∑iYijEqlnβik∑iYij+Eqlnθjk+constant,where ∑kρjk=1, and have the variationalθjk|ρ,β∼Gammaθjk|a0+ρjk∑i=1Dyij,b0+∑i=1DEqβik.

Simulation

Simulation scheme

We generated a single-cell data set from cell type–sorted bulk RNA-seq reference data from the DICE (Database of Immune Cell Expression, Expression quantitative trait loci and Epigenomics) project (Schmiedel et al, 2018) to evaluate the performance of ASAP and other methods with gold standard annotations. We considered the population of 13 cell types, including CD4+ T cells, CD8+ T cells, NK cells, B cells, and monocytes. We use a well-established data-generating framework implemented in the scDesign2 package—a simulation method based on a Gaussian copula-based sampling scheme (Sun et al, 2021). Gene–gene correlation structures are accurately captured in data simulated by a copula method, and the non-parametric nature of its density estimation step is more suitable for simulation experiments with little modelling assumption. Mathematical details can be found in the original work (Sun et al, 2021), but briefly, we simulated our data as follows.

For each cell type t, we generate a single-cell data matrix Y(t) based on a continuous version of the transformed data as follows:

1. Estimate gene-level mean μ and covariance between genes Σ for each cell type t.

2. For each cell j, we generate bootstrapped copula yj:

  • − Sample a Gaussian vector zj ∼ N(μt, Σt).

  • − Construct a stochastic version of the reference gene expression vector sj by sampling with replacement.

  • − Sort the bootstrapped gene expressions Sgj, that is, S(1)j < S(2)j < ⋯ < S(D)j.

  • − Assign each gene g’s gene expression Ygj in two steps:

    •   • Identify the ascending order r for Zgj

    •   • Take the bootstrapped S(r)j and set Ygj to this value.

3. Repeat the above steps until we sample the desired number of cells.

Next, we simulated a separate single-cell data matrix from a null model that captured the global gene expression pattern across all cell types in the bulk samples. To account for the background and cell type–invariant patterns, we generate the null data Y∼ by sampling a gene vector y∼ from a multinomial distribution with the empirical gene frequencies, ignoring cell-type identities. Depending on the ρ ∈ (0,1) value, we differently mix the cell type–specific foreground and the background signals: yj←yj(t)ρ+y∼j(1−ρ). Both foreground and background data were normalized to have the same sequencing depth.

Acknowledgements

This work was supported by the BC Cancer Foundation, NSERC Discovery, and Canada Research Chair Tier-2 (YP Park). We also acknowledge support from the UBC Four Year Fellowship (S Subedi).

Author Contributions

  • S Subedi: conceptualization, resources, data curation, software, formal analysis, funding acquisition, investigation, visualization, methodology, and writing—original draft, review, and editing.

  • TS Sumida: conceptualization, investigation, methodology, and writing—original draft.

  • YP Park: conceptualization, resources, data curation, software, formal analysis, supervision, funding acquisition, investigation, visualization, methodology, project administration, and writing—original draft, review, and editing.

Conflict of Interest Statement

The authors declare that they have no conflict of interest.

  • Received March 12, 2024.
  • Revision received July 29, 2024.
  • Accepted July 30, 2024.
Creative Commons logoCreative Commons logohttps://creativecommons.org/licenses/by/4.0/

This article is available under a Creative Commons License (Attribution 4.0 International, as described at https://creativecommons.org/licenses/by/4.0/).

References

  1. ↵
    1. Abdelaal T,
    2. de Raadt P,
    3. Lelieveldt BPF,
    4. Reinders MJT,
    5. Mahfouz A
    (2020) SCHNEL: Scalable clustering of high dimensional single-cell data. Bioinformatics 36: i849–i856. doi:10.1093/bioinformatics/btaa816
    OpenUrlCrossRef
  2. ↵
    1. Becht E,
    2. McInnes L,
    3. Healy J,
    4. Dutertre C-A,
    5. Kwok IWH,
    6. Ng LG,
    7. Ginhoux F,
    8. Newell EW
    (2018) Dimensionality reduction for visualizing single-cell data using UMAP. Nat Biotechnol 37: 38–44. doi:10.1038/nbt.4314
    OpenUrlCrossRefPubMed
  3. ↵
    1. Bingham E,
    2. Mannila H
    (2001) Random projection in dimensionality reduction: Applications to image and text data. In Proceedings of the Seventh ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’01, pp 245–250. New York, NY: Association for Computing Machinery.
  4. ↵
    1. Blei DM,
    2. Ng AY,
    3. Jordan MI
    (2003) Latent dirichlet allocation. J Mach Learn Res 3: 993–1022.
    OpenUrlCrossRef
  5. ↵
    1. Blondel VD,
    2. Guillaume J-L,
    3. Lambiotte R,
    4. Lefebvre E
    (2008) Fast unfolding of communities in large networks. J Stat Mech Theor Exp 2008: P10008. doi:10.1088/1742-5468/2008/10/p10008
    OpenUrlCrossRef
  6. ↵
    1. Saunders C,
    2. Grobelnik K,
    3. Gunn S,
    4. Shawe-Taylor J
    1. Blum A
    (2006) Random projection, margins, kernels, and feature-selection. In Subspace, Latent Structure and Feature Selection. SLSFS 2005. Lecture Notes in Computer Science, Saunders C, Grobelnik K, Gunn S, Shawe-Taylor J (eds), Vol. 3940. Berlin, Heidelberg: Springer.
  7. ↵
    1. Bravo González-Blas C,
    2. Minnoye L,
    3. Papasokrati D,
    4. Aibar S,
    5. Hulselmans G,
    6. Christiaens V,
    7. Davie K,
    8. Wouters J,
    9. Aerts S
    (2019) CisTopic: Cis-regulatory topic modeling on single-cell ATAC-seq data. Nat Methods 16: 397–400. doi:10.1038/s41592-019-0367-1
    OpenUrlCrossRefPubMed
  8. ↵
    1. Carbonetto P,
    2. Sarkar A,
    3. Wang Z,
    4. Stephens M
    (2021) Non-negative matrix factorization algorithms greatly improve topic model fits. arXiv. doi:10.48550/arXiv.2105.13440 (Preprint posted May 27, 2021).
    OpenUrlCrossRef
  9. ↵
    1. Carbonetto P,
    2. Luo K,
    3. Sarkar A,
    4. Hung A,
    5. Tayeb K,
    6. Pott S,
    7. Stephens M
    (2023) Interpreting structure in sequence count data with differential expression analysis allowing for grades of membership. BioRxiv. doi:10.1101/2023.03.03.531029 (Preprint posted June 23, 2023).
    OpenUrlAbstract/FREE Full Text
  10. ↵
    1. Chen H,
    2. Lareau C,
    3. Andreani T,
    4. Vinyard ME,
    5. Garcia SP,
    6. Clement K,
    7. Andrade-Navarro MA,
    8. Buenrostro JD,
    9. Pinello L
    (2019) Assessment of computational methods for the analysis of single-cell ATAC-seq data. Genome Biol 20: 241. doi:10.1186/s13059-019-1854-5
    OpenUrlCrossRefPubMed
  11. ↵
    1. Dasgupta S,
    2. Freund Y
    (2008) Random projection trees and low dimensional manifolds. In Proceedings of the Fortieth Annual ACM Symposium on Theory of Computing, STOC ’08, pp 537–546. New York, NY: ACM.
  12. ↵
    1. Demetci P,
    2. Santorella R,
    3. Sandstede B,
    4. Noble WS,
    5. Singh R
    (2022) SCOT: Single-cell multi-omics alignment with optimal transport. J Comput Biol 29: 3–18. doi:10.1089/cmb.2021.0446
    OpenUrlCrossRef
  13. ↵
    1. Dey KK,
    2. Hsiao CJ,
    3. Stephens M
    (2017) Visualizing the structure of RNA-seq expression data using grade of membership models. PLoS Genet 13: e1006599. doi:10.1371/journal.pgen.1006599
    OpenUrlCrossRef
  14. ↵
    1. Dhapola P,
    2. Rodhe J,
    3. Olofzon R,
    4. Bonald T,
    5. Erlandsson E,
    6. Soneji S,
    7. Karlsson G
    (2022) Scarf enables a highly memory-efficient analysis of large-scale single-cell genomics data. Nat Commun 13: 4616. doi:10.1038/s41467-022-32097-3
    OpenUrlCrossRef
  15. ↵
    1. Dieng AB,
    2. Ruiz FJR,
    3. Blei DM
    (2020) Topic modeling in embedding spaces. Trans Assoc Comput Linguistics 8: 439–453. doi:10.1162/tacl_a_00325
    OpenUrlCrossRef
  16. ↵
    1. Elyanow R,
    2. Dumitrascu B,
    3. Engelhardt BE,
    4. Raphael BJ
    (2020) netNMF-sc: Leveraging gene-gene interactions for imputation and dimensionality reduction in single-cell expression analysis. Genome Res 30: 195–204. doi:10.1101/gr.251603.119
    OpenUrlAbstract/FREE Full Text
  17. ↵
    1. Eraslan G,
    2. Drokhlyansky E,
    3. Anand S,
    4. Fiskin E,
    5. Subramanian A,
    6. Slyper M,
    7. Wang J,
    8. Van Wittenberghe N,
    9. Rouhana JM,
    10. Waldman J, et al.
    (2022) Single-nucleus cross-tissue molecular reference maps toward understanding disease gene function. Science 376: eabl4290. doi:10.1126/science.abl4290
    OpenUrlCrossRef
  18. ↵
    1. Franzén O,
    2. Gan L-M,
    3. Björkegren JLM
    (2019) PanglaoDB: A web server for exploration of mouse and human single-cell RNA sequencing data. Database 2019: baz046. doi:10.1093/database/baz046
    OpenUrlCrossRefPubMed
  19. ↵
    1. Gopalan P,
    2. Charlin L,
    3. Blei DM
    (2014) Content-based recommendations with Poisson factorization. In Proceedings of the 27th International Conference on Neural Information Processing Systems - Volume 2, NIPS’14, pp 3176–3184. Cambridge, MA: MIT Press.
  20. ↵
    1. Gopalan P,
    2. Hao W,
    3. Blei DM,
    4. Storey JD
    (2016) Scaling probabilistic models of genetic variation to millions of humans. Nat Genet 48: 1587–1590. doi:10.1038/ng.3710
    OpenUrlCrossRef
  21. ↵
    1. Granja JM,
    2. Corces MR,
    3. Pierce SE,
    4. Bagdatli ST,
    5. Choudhry H,
    6. Chang HY,
    7. Greenleaf WJ
    (2021) Author correction: ArchR is a scalable software package for integrative single-cell chromatin accessibility analysis. Nat Genet 53: 935. doi:10.1038/s41588-021-00850-x
    OpenUrlCrossRef
  22. ↵
    1. GTEx Consortium
    (2020) The GTEx consortium atlas of genetic regulatory effects across human tissues. Science 369: 1318–1330. doi:10.1126/science.aaz1776
    OpenUrlAbstract/FREE Full Text
  23. ↵
    1. Haghverdi L,
    2. Lun ATL,
    3. Morgan MD,
    4. Marioni JC
    (2018) Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors. Nat Biotechnol 36: 421–427. doi:10.1038/nbt.4091
    OpenUrlCrossRefPubMed
  24. ↵
    1. Hao Y,
    2. Hao S,
    3. Andersen-Nissen E,
    4. Mauck WM 3rd.,
    5. Zheng S,
    6. Butler A,
    7. Lee MJ,
    8. Wilk AJ,
    9. Darby C,
    10. Zager M, et al.
    (2021) Integrated analysis of multimodal single-cell data. Cell 184: 3573–3587.e29. doi:10.1016/j.cell.2021.04.048
    OpenUrlCrossRefPubMed
  25. ↵
    1. Heumos L,
    2. Schaar AC,
    3. Lance C,
    4. Litinetskaya A,
    5. Drost F,
    6. Zappia L,
    7. Lücken MD,
    8. Strobl DC,
    9. Henao J,
    10. Curion F, et al.
    (2023) Best practices for single-cell analysis across modalities. Nat Rev Genet 24: 550–572. doi:10.1038/s41576-023-00586-w
    OpenUrlCrossRef
  26. ↵
    1. Hie B,
    2. Bryson B,
    3. Berger B
    (2019a) Efficient integration of heterogeneous single-cell transcriptomes using scanorama. Nat Biotechnol 37: 685–691. doi:10.1038/s41587-019-0113-3
    OpenUrlCrossRefPubMed
  27. ↵
    1. Hie B,
    2. Cho H,
    3. DeMeo B,
    4. Bryson B,
    5. Berger B
    (2019b) Geometric sketching compactly summarizes the Single-Cell transcriptomic landscape. Cell Syst 8: 483–493.e7. doi:10.1016/j.cels.2019.05.003
    OpenUrlCrossRef
  28. ↵
    1. Hotelling H
    (1933) Analysis of a complex of statistical variables into principal components. J Educ Psychol 24: 417–441. doi:10.1037/h0071325
    OpenUrlCrossRef
  29. ↵
    1. Hristov BH,
    2. Bilmes JA,
    3. Noble WS
    (2022) Linking cells across single-cell modalities by synergistic matching of neighborhood structure. Bioinformatics 38: ii148–ii154. doi:10.1093/bioinformatics/btac481
    OpenUrlCrossRef
  30. ↵
    1. Johnson WB,
    2. Lindenstrauss J
    (1984) Extensions of lipschitz mappings into a hilbert space. Conf Mod Anal Probab: 189–206. doi:10.1090/conm/026/737400
    OpenUrlCrossRef
  31. ↵
    1. Jolliffe IT
    1. Jolliffe IT
    (1986) Principal component analysis and factor analysis. In Principal Component Analysis, Jolliffe IT (ed), pp 115–128. New York, NY: Springer New York.
  32. ↵
    1. Kazwini NE,
    2. Sanguinetti G
    (2024) SHARE-topic: Bayesian interpretable modeling of single-cell multi-omic data. Genome Biol 25: 55. doi:10.1186/s13059-024-03180-3
    OpenUrlCrossRef
  33. ↵
    1. Kharchenko PV
    (2021) The triumphs and limitations of computational methods for scRNA-seq. Nat Methods 18: 723–732. doi:10.1038/s41592-021-01171-x
    OpenUrlCrossRefPubMed
  34. ↵
    1. Korotkevich G,
    2. Sukhov V,
    3. Budin N,
    4. Shpak B,
    5. Artyomov MN,
    6. Sergushichev A
    (2021) Fast gene set enrichment analysis. BioRxiv. doi:10.1101/060012 (Preprint posted February 01, 2021).
    OpenUrlAbstract/FREE Full Text
  35. ↵
    1. Kriebel AR,
    2. Welch JD
    (2022) UINMF performs mosaic integration of single-cell multi-omic datasets using nonnegative matrix factorization. Nat Commun 13: 780. doi:10.1038/s41467-022-28431-4
    OpenUrlCrossRef
  36. ↵
    1. Levitin HM,
    2. Yuan J,
    3. Cheng YL,
    4. Ruiz F, Jr.,
    5. Bush EC,
    6. Bruce JN,
    7. Canoll P,
    8. Iavarone A,
    9. Lasorella A,
    10. Blei DM, et al.
    (2019) De novo gene signature identification from single-cell RNA-seq with hierarchical Poisson factorization. Mol Syst Biol 15: e8557. doi:10.15252/msb.20188557
    OpenUrlAbstract/FREE Full Text
  37. ↵
    1. Liu B,
    2. Li C,
    3. Li Z,
    4. Wang D,
    5. Ren X,
    6. Zhang Z
    (2020) An entropy-based metric for assessing the purity of single cell populations. Nat Commun 11: 3155. doi:10.1038/s41467-020-16904-3
    OpenUrlCrossRef
  38. ↵
    1. Moschitti A,
    2. Pang B,
    3. Daelemans W
    1. Mimno D,
    2. Lee M
    (2014) Low-dimensional embeddings for interpretable anchor-based topic inference. In Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP), Moschitti A, Pang B, Daelemans W (eds), pp 1319–1328. Doha: Association for Computational Linguistics.
  39. ↵
    1. Novembre J
    (2016) Pritchard, stephens, and donnelly on population structure. Genetics 204: 391–393. doi:10.1534/genetics.116.195164
    OpenUrlFREE Full Text
  40. ↵
    1. Pedregosa F,
    2. Varoquaux G,
    3. Gramfort A,
    4. Michel V,
    5. Thirion B,
    6. Grisel O,
    7. Blondel M,
    8. Prettenhofer P,
    9. Weiss R,
    10. Dubourg V, et al.
    (2011) Scikit-learn: Machine learning in Python. J Mach Learn Res 12: 2825–2830.
    OpenUrlCrossRefPubMed
  41. ↵
    1. Persad S,
    2. Choo Z-N,
    3. Dien C,
    4. Masilionis I,
    5. Chaligné R,
    6. Nawy T,
    7. Brown CC,
    8. Pe’er I,
    9. Setty M,
    10. Pe’er D
    (2022) SEACells: Inference of transcriptional and epigenomic cellular states from single-cell genomics data. BioRxiv. doi:10.1101/2022.04.02.486748 (Preprint posted April 03, 2022).
    OpenUrlAbstract/FREE Full Text
  42. ↵
    1. Polański K,
    2. Young MD,
    3. Miao Z,
    4. Meyer KB,
    5. Teichmann SA,
    6. Park JE
    (2020) BBKNN: Fast batch alignment of single cell transcriptomes. Bioinformatics 36: 964–965. doi:10.1093/bioinformatics/btz625
    OpenUrlCrossRefPubMed
  43. ↵
    1. Pritchard JK,
    2. Stephens M,
    3. Donnelly P
    (2000) Inference of population structure using multilocus Genotype data. Genetics 155: 945–959. doi:10.1093/genetics/155.2.945
    OpenUrlAbstract/FREE Full Text
  44. ↵
    1. Rand WM
    (1971) Objective criteria for the evaluation of clustering methods. J Am Stat Assoc 66: 846–850. doi:10.2307/2284239
    OpenUrlCrossRef
  45. ↵
    1. Schmiedel BJ,
    2. Singh D,
    3. Madrigal A,
    4. Valdovino-Gonzalez AG,
    5. White BM,
    6. Zapardiel-Gonzalo J,
    7. Ha B,
    8. Altay G,
    9. Greenbaum JA,
    10. McVicker G, et al.
    (2018) Impact of genetic polymorphisms on human immune cell gene expression. Cell 175: 1701–1715.e16. doi:10.1016/j.cell.2018.10.022
    OpenUrlCrossRefPubMed
  46. ↵
    1. Stoeckius M,
    2. Hafemeister C,
    3. Stephenson W,
    4. Houck-Loomis B,
    5. Chattopadhyay PK,
    6. Swerdlow H,
    7. Satija R,
    8. Smibert P
    (2017) Simultaneous epitope and transcriptome measurement in single cells. Nat Methods 14: 865–868. doi:10.1038/nmeth.4380
    OpenUrlCrossRefPubMed
  47. ↵
    1. Stuart T,
    2. Satija R
    (2019) Integrative single-cell analysis. Nat Rev Genet 20: 257–272. doi:10.1038/s41576-019-0093-7
    OpenUrlCrossRefPubMed
  48. ↵
    1. Subedi S,
    2. Park YP
    (2023) Single-cell pair-wise relationships untangled by composite embedding model. iScience 26: 106025. doi:10.1016/j.isci.2023.106025
    OpenUrlCrossRef
  49. ↵
    1. Subramanian A,
    2. Tamayo P,
    3. Mootha VK,
    4. Mukherjee S,
    5. Ebert BL,
    6. Gillette MA,
    7. Paulovich A,
    8. Pomeroy SL,
    9. Golub TR,
    10. Lander ES, et al.
    (2005) Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A 102: 15545–15550. doi:10.1073/pnas.0506580102
    OpenUrlAbstract/FREE Full Text
  50. ↵
    1. Sun S,
    2. Zhu J,
    3. Ma Y,
    4. Zhou X
    (2019) Accuracy, robustness and scalability of dimensionality reduction methods for single-cell RNA-seq analysis. Genome Biol 20: 269. doi:10.1186/s13059-019-1898-6
    OpenUrlCrossRef
  51. ↵
    1. Sun T,
    2. Song D,
    3. Li WV,
    4. Li JJ
    (2021) ScDesign2: A transparent simulator that generates high-fidelity single-cell gene expression count data with gene correlations captured. Genome Biol 22: 163. doi:10.1186/s13059-021-02367-2
    OpenUrlCrossRef
  52. ↵
    1. Tabula Sapiens Consortium,
    2. Jones RC,
    3. Karkanias J,
    4. Krasnow MA,
    5. Pisco AO,
    6. Quake SR,
    7. Salzman J,
    8. Yosef N,
    9. Bulthaup B,
    10. Brown P, et al.
    (2022) The tabula sapiens: A multiple-organ, single-cell transcriptomic atlas of humans. Science 376: eabl4896. doi:10.1126/science.abl4896
    OpenUrlCrossRefPubMed
  53. ↵
    1. Townes FW,
    2. Engelhardt BE
    (2023) Nonnegative spatial factorization applied to spatial genomics. Nat Methods 20: 229–238. doi:10.1038/s41592-022-01687-w
    OpenUrlCrossRef
  54. ↵
    1. Traag V,
    2. Waltman L,
    3. van Eck NJ
    (2018) From louvain to leiden: Guaranteeing well-connected communities. arXiv. doi:10.48550/arXiv.1810.08473 (Preprint posted October 19, 2018).
    OpenUrlCrossRef
  55. ↵
    1. Tsuyuzaki K,
    2. Sato H,
    3. Sato K,
    4. Nikaido I
    (2020) Benchmarking principal component analysis for large-scale single-cell RNA-sequencing. Genome Biol 21: 9. doi:10.1186/s13059-019-1900-3
    OpenUrlCrossRef
  56. ↵
    1. Vinh NX,
    2. Epps J,
    3. Bailey J
    (2010) Information theoretic measures for clusterings comparison: Variants, properties, normalization and correction for chance. J Mach Learn Res 11: 2837–2854.
    OpenUrlCrossRef
  57. ↵
    1. Waltman L,
    2. van Eck NJ
    (2013) A smart local moving algorithm for large-scale modularity-based community detection. Eur Phys J B 86: 471. doi:10.1140/epjb/e2013-40829-0
    OpenUrlCrossRef
  58. ↵
    1. Wan S,
    2. Kim J,
    3. Won KJ
    (2020) SHARP: Hyperfast and accurate processing of single-cell RNA-seq data via ensemble random projection. Genome Res 30: 205–213. doi:10.1101/gr.254557.119
    OpenUrlAbstract/FREE Full Text
  59. ↵
    1. Wei N,
    2. Nie Y,
    3. Liu L,
    4. Zheng X,
    5. Wu H-J
    (2022) Secuer: Ultrafast, scalable and accurate clustering of single-cell RNA-seq data. PLoS Comput Biol 18: e1010753. doi:10.1371/journal.pcbi.1010753
    OpenUrlCrossRef
  60. ↵
    1. Welch JD,
    2. Kozareva V,
    3. Ferreira A,
    4. Vanderburg C,
    5. Martin C,
    6. Macosko EZ
    (2019) Single-cell multi-omic integration compares and contrasts features of brain cell identity. Cell 177: 1873–1887.e17. doi:10.1016/j.cell.2019.05.006
    OpenUrlCrossRefPubMed
  61. ↵
    1. Wolf FA,
    2. Angerer P,
    3. Theis FJ
    (2018) SCANPY: Large-scale single-cell gene expression data analysis. Genome Biol 19: 15. doi:10.1186/s13059-017-1382-0
    OpenUrlCrossRefPubMed
  62. ↵
    1. Wu SZ,
    2. Al-Eryani G,
    3. Roden DL,
    4. Junankar S,
    5. Harvey K,
    6. Andersson A,
    7. Thennavan A,
    8. Wang C,
    9. Torpy JR,
    10. Bartonicek N, et al.
    (2021) A single-cell and spatially resolved atlas of human breast cancers. Nat Genet 53: 1334–1347. doi:10.1038/s41588-021-00911-1
    OpenUrlCrossRefPubMed
  63. ↵
    1. Zhang Y,
    2. Khalilitousi M,
    3. Park YP
    (2023) Unraveling dynamically encoded latent transcriptomic patterns in pancreatic cancer cells by topic modeling. Cell Genomics 3: 100388. doi:10.1016/j.xgen.2023.100388
    OpenUrlCrossRef
  64. ↵
    1. Zhang K,
    2. Zemke NR,
    3. Armand EJ,
    4. Ren B
    (2024) A fast, scalable and versatile tool for analysis of single-cell omics data. Nat Methods 21: 217–227. doi:10.1038/s41592-023-02139-9
    OpenUrlCrossRef
  65. ↵
    1. Zhao X,
    2. Li H,
    3. Lyu S,
    4. Zhai J,
    5. Ji Z,
    6. Zhang Z,
    7. Zhang X,
    8. Liu Z,
    9. Wang H,
    10. Xu J, et al.
    (2021a) Single-cell transcriptomics reveals heterogeneous progression and EGFR activation in pancreatic adenosquamous carcinoma. Int J Biol Sci 17: 2590–2605. doi:10.7150/ijbs.58886
    OpenUrlCrossRef
  66. ↵
    1. Zhao Y,
    2. Cai H,
    3. Zhang Z,
    4. Tang J,
    5. Li Y
    (2021b) Learning interpretable cellular and gene signature embeddings from single-cell transcriptomic data. Nat Commun 12: 5261. doi:10.1038/s41467-021-25534-2
    OpenUrlCrossRef
Back to top
Download PDF
Email Article
Citation Tools
Single-cell topic analysis by ASAP
Sishir Subedi, Tomokazu S Sumida, Yongjin P Park
Life Science Alliance Aug 2024, 7 (10) e202402713; DOI: 10.26508/lsa.202402713
Twitter logo Facebook logo Mendeley logo
Issue Cover

In this Issue

Volume 7, No. 10
October 2024
Advertisement

Jump to section

Subjects

Related Articles

  • No related articles found.

Cited By...

  • No citing articles found.

More in this TOC Section

Similar Articles