Author response: Self-assembling manifolds in single-cell RNA sequencing data (2019)

Article Figures and data Abstract eLife digest Introduction Results Discussion Materials and methods Data availability References Decision letter Author response Article and author information Metrics Abstract Single-cell RNA sequencing has spurred the development of computational methods that enable researchers to classify cell types, delineate developmental trajectories, and measure molecular responses to external perturbations. Many of these technologies rely on their ability to detect genes whose cell-to-cell variations arise from the biological processes of interest rather than transcriptional or technical noise. However, for datasets in which the biologically relevant differences between cells are subtle, identifying these genes is challenging. We present the self-assembling manifold (SAM) algorithm, an iterative soft feature selection strategy to quantify gene relevance and improve dimensionality reduction. We demonstrate its advantages over other state-of-the-art methods with experimental validation in identifying novel stem cell populations of Schistosoma mansoni, a prevalent parasite that infects hundreds of millions of people. Extending our analysis to a total of 56 datasets, we show that SAM is generalizable and consistently outperforms other methods in a variety of biological and quantitative benchmarks. https://doi.org/10.7554/eLife.48994.001 eLife digest New technologies have enabled scientists to closely examine the activity of individual cells. One increasingly popular technique to do this is called single-cell RNA-sequencing and it relies on the fact that although all cells in an organism carry the same DNA, different cell types use different genes. This technique is powerful but can struggle to identify meaningful distinctions between cell types, especially when the differences are subtle. In single-cell RNA-sequencing, the messenger RNA (mRNA) copied from each gene is collected and counted, and usually the more a gene is copied the more active it is. Differences in gene activity (also called gene expression) between two cells often imply that they are different types of cells. However, since only an infinitesimal amount of mRNAs can be collected from a single cell, the counting is often inaccurate. In addition, the transient changes in gene expression can make cells of the same type have different gene expressions. These factors make it challenging to determine what genes are informative for distinguishing between cell types. To address this problem, Tarashansky et al. have developed a computational approach called Self-Assembling Manifold (or SAM for short) to identify differences in gene expression that can lead to a better classification of cell types. First, SAM groups the cells randomly and looks for genes with different expression patterns between the groups. By looking at differences between groups instead of differences between individual cells, SAM is ‘averaging out’ individual differences within groups. SAM then uses the resulting information to re-classify the cells and start the process over again, taking the new groups and finding differences between them. SAM repeats these steps until the classification stops changing and becomes stable. SAM does not require any existing knowledge about cell types or gene expression, meaning it is unbiased and widely applicable. To test the usefulness of the algorithm, Tarashansky et al. used SAM to identify new cell types in the medically important parasitic worm Schistosoma mansoni, which infects hundreds of millions of people worldwide every year. SAM can tell cell types apart better than existing approaches, and it can find meaningful differences in systems with a lot of meaningless variability as demonstrated by evaluating SAM’s performance on 55 other datasets. The potential applications of this approach are many, including the creation of detailed cell atlases recording the different types of cells throughout entire organisms. https://doi.org/10.7554/eLife.48994.002 Introduction Single-cell RNA sequencing (scRNAseq) datasets typically contain tens of thousands of genes, although many of them may not be informative for differentiating between cell types or states. Feature selection is thus commonly used to select a subset of genes prior to downstream analyses, such as manifold reconstruction and cell clustering (Crow et al., 2018; Satija et al., 2015; Vallejos et al., 2015). However, current approaches have two major limitations. First, feature selection methods filter genes based on arbitrarily or empirically chosen thresholds, small changes in which may result in different gene sets (Vallejos et al., 2017). In addition, the selection of features typically operates under the assumption that genes with highly variable expression between individual cells capture biologically meaningful variation. Because single-cell transcriptomes are inevitably contaminated by a combination of random transcriptional and technical noise (Grün et al., 2014), the variation in biologically relevant genes may be hard to distinguish from the background noise, especially when the differences between cell populations are subtle. Resolving these differences, or ‘signals’, is essential to study a variety of biological problems, including identifying cell subtypes (Olsson et al., 2016; Treutlein et al., 2014; Lönnberg et al., 2017; Fincher et al., 2018; Baron et al., 2016; Schwalie et al., 2018) and quantifying the effects of molecular perturbations to otherwise homogeneous populations of cells (Lane et al., 2017). In such datasets, only a small fraction of the genes, and therefore only a small fraction of the total variation, may contain the signals relevant for distinguishing cell types or cell states. Choosing these features without a priori knowledge remains an unmet computational challenge. The second limitation is that existing methods have been almost exclusively benchmarked on well-annotated, gold standard datasets with clearly distinguishable cell types (Wang et al., 2017; Kiselev et al., 2017; Duò et al., 2019; Bahlo et al., 2018). These datasets are not informative for distinguishing the performance between methods, because the differences between cell types are relatively straightforward to detect. However, evaluating the performance of feature selection and/or dimensionality reduction methods on datasets with more subtle signals is difficult as their ground truth labels are typically ambiguous or nonexistent. To overcome the shortcomings of current feature selection approaches, here, we introduce the Self-Assembling Manifold (SAM) method, an unsupervised, ‘soft feature selection’ algorithm that iteratively rescales gene expressions to refine a nearest neighbor graph of cells until the graph converges to a stable solution. At each iteration, SAM assigns more weight to genes that are spatially variable across the constructed graph, and this weighted gene expression is then used to improve the next nearest neighbor assignment. SAM presents two advantages: it rescales all genes according to their weights, solving the problem of thresholding, and it prioritizes genes that are variable across the intrinsic manifold of the data rather than selecting genes that are variable across individual cells. In order to better distinguish the performance between methods, we define a network sensitivity measure to identify datasets with subtle signals. With limited annotations in most high-sensitivity datasets, we introduce unsupervised graph-based metrics to quantify the degree of structure within the reconstructed manifolds for comparison between methods. In addition, we perform benchmarking using known ground truth labels on simulated datasets spanning a wide range of sensitivities by introducing increasing levels of noise to well-annotated datasets. These analyses reveal that SAM consistently improves feature selection and cell clustering. To demonstrate the utility of SAM in practice, we provide an in-depth analysis of two datasets that are challenging to analyze using existing methods: stem cells in a human parasitic worm, Schistosoma, and activated macrophages (Lane et al., 2017). We show that SAM can capture novel biology undetectable by other methods and validate these results with experimental evidence. Results The SAM algorithm The SAM algorithm begins with a random k-nearest neighbor (kNN) graph and averages the expression of each cell with its k nearest neighbors: C=1kNE, where N is the directed adjacency matrix and E is the gene expression matrix (Figure 1a). For each gene i, SAM computes a spatial dispersion factor of the averaged expressions Ci, which measures variation across neighborhoods of cells rather than individual cells (Materials and methods). These dispersions are used to calculate the gene weights, which then rescale the expression matrix: E^=EWD, where WD is a diagonal matrix with gene weights along the diagonal. Using the rescaled expressions E^, we compute a pairwise cell distance matrix and update the assignment of each cell’s k-nearest neighbors accordingly. This cycle continues until the gene weights converge. To demonstrate the implementation and utility of SAM, below we analyze a challenging dataset comprised of a few hundred relatively homogeneous stem cells isolated from Schistosoma mansoni (Figure 1—figure supplement 1), a widespread human pathogen (Hoffmann et al., 2014). Using a protocol we have established previously (Wang et al., 2018), these cells were collected by sorting dividing cells from juvenile parasites harvested from their mouse hosts at 2.5 weeks post infection. At this stage, the parasites use an abundant stem cell population (~15–20% of the total number of cells) for rapid organogenesis and growth (Wang et al., 2013; Wang et al., 2018). Testing several existing methods (Wang et al., 2017; Kiselev et al., 2017; Satija et al., 2015), we found that they were not able to identify distinct cell populations in this dataset. In contrast, SAM finds a stable solution independent of initial conditions (Figure 1b). A graph structure with clearly separated cell populations self-assembles through the iterative process (Figure 1c). In parallel, the gene weights converge onto the final weight vector. Eventually, only a small fraction of genes (~1%) are strongly weighted and useful for separating cell clusters, reflecting the inherent difficulty of analyzing this dataset. Figure 1 with 2 supplements see all Download asset Open asset The SAM algorithm. (a) SAM starts with a randomly initialized kNN adjacency matrix and iterates to refine the adjacency matrix and gene weight vector until convergence. (b) Root mean square error (RMSE) of the gene weights (top) and the fraction of different edges of the nearest-neighbor adjacency matrices (bottom) between adjacent iterations (blue) and between independent runs at the same iteration (orange) to show that SAM converges to the same solution regardless of initial conditions. The differences between the gene weights and nearest-neighbor graphs from independent runs are relatively small, indicating that SAM converges to the same solution through similar paths. (c) Graph structures and gene weights of the schistosome stem cell data converging to the final output over the course of 10 iterations (i denotes iteration number). Top: nodes are cells and edges connect neighbors. Nodes are color-coded according to the final clusters. Bottom: weights are sorted according to the final gene rankings. (d) Network properties iteratively improve for the graphs reconstructed from the original data (red) but not on the randomly shuffled data (blue). The network properties converge to the same values when initializing SAM with the Seurat-reconstructed graph instead of a random graph (yellow). Dashed lines: metrics measured from the Seurat-reconstructed graphs. https://doi.org/10.7554/eLife.48994.003 Figure 1d shows that SAM iteratively improves a series of graph characteristics, including the network-average clustering coefficient (NACC), modularity, and Euclidean norm of the spatial dispersions (Materials and methods). The NACC and modularity quantify the degree of structure within the graphs – graphs with high NACC and modularity have regions of high density separated by regions of low density. The dispersion quantifies the spatial organization of gene expression – the higher the spatial dispersion the less uniformly distributed the gene expressions are along the graph. The final graph metrics are independent of initial conditions, which can start from a random graph or the output of an existing manifold reconstruction algorithm (e.g. Seurat, Satija et al., 2015). Importantly, we verified that SAM does not artificially boost these metrics in data that lack inherent structure: when applying SAM to a randomly shuffled expression matrix, none of these metrics increased from the random initial conditions. SAM identifies novel subpopulations within schistosome stem cells Visualizing the converged graph in two dimensions using Uniform Manifold Approximation and Projection (Becht et al., 2019), we find that cells can be separated into three main populations, with Louvain clustering (Blondel et al., 2008) further splitting one of these clusters into two subpopulations (Figure 2a). In contrast, other commonly used dimensionality reduction methods, such as principal component analysis (PCA), Seurat (Satija et al., 2015), and SIMLR (Wang et al., 2017), failed to distinguish these cell populations (see Materials and methods for the selection of algorithms for comparison). Clustering the Seurat graph using Louvain clustering still results in a low-modularity partition and poor correspondence to the SAM cluster assignments. Figure 2 with 1 supplement see all Download asset Open asset SAM identifies novel subpopulations within schistosome stem cells. (a) UMAP projections of the manifolds reconstructed by SAM, PCA, and Seurat. SIMLR outputs its own 2D projection based on its constructed similarity matrix using a modified version of t-SNE. The schistosome cells are color-coded by the stem cell subpopulations μ, δ’, εɑ, and εβ determined by Louvain clustering. (b) UMAP projections with gene expressions of subpopulation-specific markers (eledh, nanos-2, cabp, astf, bhlh,) and a ubiquitous stem cell marker, ago2-1, overlaid. Insets: magnified views of the expressing populations. (c) FISH of cabp and EdU labeling of dividing stem cells in juvenile parasites at 2.5 weeks post-infection show that μ-cells (cabp+EdU+, arrowheads) are close to the parasite surface and beneath a layer of post-mitotic cabp+ cells. Dashed outline: parasite surface. Right: magnified views of the boxed region. (d) FISH of cabp and a set of canonical muscle markers, troponin, myosin, tropomyosin, and collagen, shows colocalization in post-mitotic cabp+ cells. Images in (c–d) are single confocal slices. (e) FISH of astf and bhlh shows their orthogonal expression in adjacent EdU+ cells (arrowheads). Bottom: magnified views of the boxed region. Image is a maximum intensity projection of a confocal stack with a thickness of 12 µm. (f) UMAP projection of stem cells isolated from juveniles at 2.5 and 3.5 weeks post-infection. Cell subpopulation assignments based on marker gene expressions are specified. Right: a magnified view to show the mapping of εɑ- and εβ-cells. (g) Standardized dispersions as calculated by Seurat plotted vs. the SAM gene weights. (h) SC3 AUROC scores plotted vs. the SAM gene weights. Error bars indicate the standard deviation of SC3 AUROC scores between trials using different chosen numbers of clusters. In (g) and (h), the top 20 genes specific to each subpopulation are colored according to the color scheme used in (a). https://doi.org/10.7554/eLife.48994.006 Supplementary file 1 lists genes with high SAM weights, which includes most markers previously implicated to be enriched in subsets of schistosome stem cells (Wang et al., 2013; Wang et al., 2018). Figure 2b shows that the three populations include previously characterized δ′-cells, which specifically express an RNA binding protein nanos-2 (Smp_051920), and ε-cells, which are marked by the expression of eledh (eled, Smp_041540) (Wang et al., 2018). More importantly, SAM reveals a novel stem cell population, μ, comprising ~30% of all sequenced cells (μ denotes muscle progenitors as discussed below). While μ-cells express ubiquitous stem cells markers (e.g. ago2-1, Smp_179320; cyclin B, Smp_082490) and cell cycle regulators (Figure 2—figure supplement 1a) (Collins et al., 2013; Wang et al., 2013; Wang et al., 2018), they are also strongly enriched for a large set of genes, with a calcium binding protein (cabp, Smp_005350), an actin protein (Smp_161920), an annexin homolog (Smp_074140), a helix-loop-helix transcription factor (dhand, Smp_062490), and a phosphatase (dusp10, Smp_034500) as the most specific markers of this population in comparison to other stem cells (Figure 2—figure supplement 1b). Fluorescent in-situ hybridization (FISH) in conjunction with EdU labeling of dividing cells reveals that μ-cells (cabp+EdU+) are distributed near the parasite surface right beneath a layer of post-mitotic differentiated cells that also express cabp (Figure 2c). Close to the parasite surface, there are two major cell types intertwined in space: the skin-like tegumental cells and the body wall muscle cells. However, μ-cells express none of the recently identified markers in tegumental progenitors (Wendt et al., 2018), suggesting that they may be associated with the muscle lineage. To test this idea, we performed double FISH experiments and observed in post-mitotic cabp+ cells the coexpression of a set of canonical muscle markers (Witchley et al., 2013), including tropomyosin (Smp_031770), myosin (Smp_045220), troponin (Smp_018250), and collagen (Smp_170340) (Figure 2d). These results suggest that cabp may mark the parasite body wall muscles and μ-cells are likely muscle progenitors, although functional validation is required to support this observation. Why the juvenile parasites maintain such an active pool of muscle progenitors will be an important question for future studies. In addition, SAM identifies two subpopulations among ε-cells: εɑ-cells that are highly enriched for an aschaete-scute transcription factor (astf, Smp_142120), and εβ-cells that abundantly express another basic helix-loop-helix protein (bhlh, Smp_087310) (Figure 2b, right panels). FISH experiments confirm these cells to be in close spatial proximity but with no coexpression of astf and bhlh (Figure 2e). Moreover, we observed with FISH that there are fewer astf+ cells in larger, more matured juveniles, suggesting εɑ-cells are a dynamic population during development. To verify this observation, we sequenced another ~370 stem cell from juveniles at a later developmental time point (3.5 weeks post infection). After correcting for batch effects in the combined 2.5- and 3.5- week datasets using the mutual nearest neighbors (MNN) algorithm (Haghverdi et al., 2018), we find that δ′-, μ-, and εβ-cells remain relatively constant throughout both time points, whereas εɑ-cells comprise a significantly smaller fraction of the stem cells at 3.5 weeks (7%) compared to 21% at 2.5 weeks (Figure 2f). Taken together, these analyses demonstrate that SAM can identify experimentally validated stem cell populations that are previously subtle to using other methods but are closely associated with the schistosome development. The between SAM and other methods in they select genes for manifold SAM prioritizes genes with variable expressions across neighborhoods of cells rather than individual cells as in other methods (e.g. Figure shows that genes with high dispersion across individual cells often have low SAM weights, indicating that these highly variable genes are to the between cells. methods (e.g. Kiselev et al., identify marker genes based on gene expression between cell clusters, but this approach when cell cluster assignment is especially when cell groups are difficult to or SC3 failed in the as it there to be only one cluster in the schistosome dataset. After we increased the number of clusters, SC3 a few of the marker genes associated with only one in Figure of the populations by changing the number of clusters in different and large variability in SC3 scores for its top genes. SAM outperforms other state-of-the-art methods in quantitative benchmarking we the of SAM by benchmarking its performance state-of-the-art analysis methods on a large of datasets. We on three methods, that Seurat, and as they are unsupervised, have been and were to other methods through benchmarking et al., 2017; Wang et al., 2017; Duò et al., 2019; Bahlo et al., 2018; et al., The to select algorithms for comparison are in Materials and methods. We datasets file that have annotations to the of SAM in cell clusters. of these datasets are of cells, as their subpopulations have been characterized with known marker genes et al., For of the datasets, SAM has the a measure of clustering and with to the annotations (Figure the Baron datasets, SAM and Seurat perform with near clustering whereas SC3 and SIMLR to and the number of clusters, Supplementary file lists the clustering scores for each and for each cell type in the benchmarking datasets (Materials and methods). SC3 and SIMLR struggle to cluster cell types that large of the such as the and cells in the datasets. While Seurat on the Baron datasets, it to identify cells in the Wang and datasets when with although its performance is to its clustering (Materials and methods). We that this is to perform on an experimental dataset with no ground truth with Seurat has than or to that of SAM on all datasets. Figure with 1 supplement see all Download asset Open asset SAM improves clustering and (a) of cluster assignment by on datasets Right: differences between the number of clusters found by each and the number of clusters differences indicate more clustering. denotes Seurat analysis using that (b) of gene weights output by SAM averaged across runs with random initial conditions for 56 datasets (blue) and simulated datasets with no intrinsic structure Materials and methods). (c) of SAM, and Seurat as a of the number of cells in each dataset. SC3 and SIMLR were not on datasets with cells as the time 20 SAM converges to the same set of gene weights for all datasets (Figure Figure supplement 1a) and its performance is to the of and random initial conditions (Figure supplement In contrast, applying SAM to randomly datasets (Materials and the resulting gene weights are highly across random initial conditions (Figure indicating that SAM does not converge to a stable solution on datasets with no intrinsic the of SAM is similar to that of Seurat, of analyzing hundreds of thousands of cells in (Figure whereas SIMLR and SC3 are of and thus from further benchmarking which the analysis of many more datasets. Because the benchmarking datasets are all comprised of clearly distinguishable cell types, they may not the performance of methods on other datasets that contain cell populations that are only To identify such datasets, we introduce a network sensitivity that quantifies the changes in the cell-to-cell when randomly selecting a subset of features from the gene expression matrices (Materials and methods). network sensitivity that changes to the features strongly the resulting that are to the features to datasets that have many signals or genes the network In the datasets we file all used benchmarking datasets have sensitivities whereas the schistosome which we have to be challenging to has the sensitivity (Figure The fraction of genes with large SAM weights is with the network suggesting that the biologically relevant variation in datasets with high sensitivity is by relatively fewer genes (Figure all 56 datasets, we found that SAM improves the modularity, and spatial organization of gene expression across the graph in comparison to Seurat as the datasets increasingly (Figure Figure Download asset Open asset SAM improves the analysis of datasets with network (a) Network sensitivity of all 56 datasets in the benchmarking datasets used in Figure measures the of a dataset to changes in which features are (Materials and methods). (b) The network sensitivity plotted the fraction of genes with SAM weight than with coefficient in the (c) of SAM over Seurat for modularity, and spatial dispersion with to sensitivity for all 56 datasets. These are with network sensitivity with in the of each the clustering for the highly datasets, is because many of them have or cell type we use the well-annotated benchmarking datasets to data across a wide of For we the data by randomly increasing of the gene expressions. by the dataset et al., 2015), Figure shows that the sensitivity along with the SAM’s scores only as the whereas performance when with A similar observed between SAM and Seurat with the modularity, and spatial Importantly, the genes with high SAM weights into Seurat its performance across all indicating that SAM is able to consistently capture the genes relevant to the structure of the data with increasing levels of noise and the of its feature selection strategy compared to the approach used by Seurat. These to all benchmarking datasets, by the under the of the metrics with to (Figure Figure Download asset Open asset feature selection improves cell clustering and manifold (a) Network modularity, and spatial dispersion with to of the in which we randomly of the data from to of the total number of (Materials and methods). is compared between SAM Seurat Seurat with and Seurat with the SAM genes Error bars indicate the standard across 10 The for with no bars are small to be (b) of the under of the metrics in (a) with to data for all datasets. Error bars indicate the standard across 10 The for data with no error bars are small to be SAM clusters macrophages by their with We next another dataset to show that SAM can biologically meaningful information that other methods