Bonsai reconstructs tree representations for distortion-free visualization and exploration of high-dimensional data

Single-cell omics methods promise to transform our understanding of the gene regulatory processes underlying cellular dynamics. Major efforts have measured the states of single cells using scRNA-seq for mRNA expression levels¹ and scATAC-seq for chromatin accessibility^2,3, including atlases of model organisms^4,5,6.

However, a consensus on how to analyze such data is yet to emerge; more than 1,750 scRNA-seq analysis tools were published in the last 8 years⁷, almost as many as the total number of scRNA-seq papers (Supplementary Fig. 1)^7,8. Analysis of these data not only faces the major technical challenges that single-cell omics data are very high-dimensional and sparse, in addition to having highly heterogeneous noise levels that vary over several orders of magnitude across measurements; the more fundamental problem is that we do not know what structure to expect in the data. For example, we do not know to what extent cells occur in discrete ‘cell types’ versus being distributed on a continuous manifold or what topology and dimensionality this manifold of cell states may have. Distributions of points in high-dimensional spaces can have extremely complex structures for which we lack intuition.

There is, thus, an urgent need for methods that can explore the complex high-dimensional structure of these data. This requires methods that can visualize local and global relationships in the data without distortion and without hallucinating structure that does not exist. Unfortunately, currently popular methods such as t-distributed stochastic neighbor embedding (t-SNE)⁹ and uniform manifold approximation and projection (UMAP)¹⁰ are well known to fail in this regard^11,12,13. In the absence of methods that reliably represent the structure in the data without distortion, the current practice in the field is to implement visualization methods in an almost trial-and-error manner^14,15,16, tweaking their tunable parameters until the visualizations match prior biological knowledge or preconceived expectations. However, this practice hinders making truly novel observations or falsifying strongly held beliefs.

Here we show that all these challenges can be overcome by representing single-cell omics data using tree structures. First, there are many examples where relationships between high-dimensional objects were efficiently captured by hierarchical representations (for example, the early taxonomies of organisms¹⁷). Trees are used ubiquitously to describe relationships between DNA sequences¹⁸ and hierarchical structures have also proven effective for representing high-dimensional objects in machine learning^19,20.

Second, as cells from a single organism are related through a lineage tree of cell divisions, gene expression patterns of single cells have in fact diverged along the branches of a tree structure. In addition, we show that, because of a feature of high-dimensional spaces that we call the ‘blessing of dimensionality’, distances between objects in high-dimensional spaces can generically be accurately represented along the branches of a tree. As trees can always be displayed in two dimensions, this allows for distortion-free visualizations of the relationships between cells.

Although conceptually similar to hierarchical clustering²¹ and the phylogeny reconstruction problem²², developing a method for relating a set of single-cell gene expression states into a most likely tree structure faced several challenges. These included developing a probability model for movement through the high-dimensional gene expression space, using this to derive a likelihood of a given tree given the expression states of cells at its leaves and developing methods for rigorously accounting for the complex heterogeneous noise properties of scRNA-seq data. In addition, we had to develop algorithms that efficiently find the maximum-likelihood tree for datasets with many thousands of cells.

We here present Bonsai, a Bayesian method that takes any set of objects with estimated coordinates in a high-dimensional continuous space, together with individual error bars on each estimated coordinate of each object, and reconstructs the most likely tree structure relating the objects at its leaves. Bonsai is derived from first principles with minimal assumptions and without any tunable parameters. Using extensive tests on real and simulated scRNA-seq data, we show that, in contrast to existing methods, Bonsai accurately represents the structure in the data on all scales. Rather than just clustering cells into different ‘cell types’, Bonsai not only infers the trajectories along which gene expression states have differentiated but also accurately preserves pairwise distances between all cells at all scales. Moreover, applying Bonsai to scRNA-seq data of cord-blood cells, it not only recovers the known differentiation hierarchy of blood cell types but also uncovers lineage relationships that, as far as we know, are previously undescribed.

In addition to providing Bonsai as a standalone tool²³, we implemented an automated pipeline for scRNA-seq analysis as a webserver (https://bonsai.unibas.ch) that starts from an mRNA count matrix, normalizes the data using Sanity²⁴, identifies clusters of statistically indistinguishable cells using Cellstates²⁵ and then runs Bonsai to reconstruct a tree relating the cells. Lastly, to facilitate exploring Bonsai’s results, we also provide Bonsai-scout, an interactive app that can be used to view the tree using different layouts, to zoom in on specific parts, to define clusters of cells by finding maximally separated clades, to overlay gene expression on the tree and to find marker genes that most distinguish different clades of cells. Example results of our integrated scRNA-seq pipeline on cell atlases of human and mouse^4,5 are provided as community resources.

As it is generally impossible to represent the distances between objects in a high-dimensional space by distances in a two-dimensional image, currently popular visualization methods focus solely on preserving nearest-neighbor relationships at a particular chosen scale, for example, such that the k-nearest-neighbor cells in the high-dimensional space are also near in the two-dimensional image²⁶. In contrast, we here aim to provide visualizations that capture the full structure in the high-dimensional data by reconstructing a tree with cells at its leaves such that, for each pair of cells, the path along the branches of the tree corresponds to the most likely trajectory through gene expression space connecting the cells. Equivalently, going out from the root of the tree to the cells at the leaves, the trajectories along the branches capture the most likely trajectories along which the gene expression states have diverged into the states at the leaves.

To define the likelihood of a tree, we need to assign probabilities to the gene expression changes along its branches and this may seem like a poorly defined problem. It is often imagined that gene expression trajectories are constrained to some manifold embedded in the high-dimensional gene expression space but we know virtually nothing about the structure or dimensionality of this manifold. However, as any smooth manifold is locally isomorphic to flat Euclidean space, at least for nearby cells, their Euclidean distance in gene expression space may be a good approximation to their distance along the manifold. As we outline in the next section, if we assume that cells have been sampled sufficiently densely that, for each individual branch of the tree, the Euclidean distance in gene expression space is a good approximation of the distance along the unknown underlying manifold, this suffices to uniquely define the likelihood of any tree connecting the gene expression states of the cells at the leaves.

Instead of defining the structure in the data in terms of the most likely trajectories through gene expression space connecting the cells, it is more common to define the quality of a visualization by how well it conserves the true distances between the cells. This raises the following question: To what extent does the tree of most likely trajectories also preserve cell-to-cell distances? Strikingly, we see below that, for high-dimensional data, the reconstructed trees in fact also accurately conserve cell-to-cell distances. We hypothesize that this is a consequence of the fact that, in high-dimensional spaces, virtually all directions are mutually orthogonal. If the underlying manifold along which gene expression states diverge is sufficiently high-dimensional, then the movements along the branches of the tree are all mutually orthogonal so that pairwise squared Euclidean distances sum along the tree’s branches. Note that this is reminiscent of what happens in the evolution of DNA sequences. Even though the evolution of protein-coding sequences is constrained to the complex subspace of sequences that code for a functional protein, this space is sufficiently high-dimensional that Hamming distances between pairs of sequences generally accurately reflect the time since the sequences diverged from a common ancestor. More surprisingly, we see below that this preservation of distances even extends to high-dimensional datasets that were not in fact generated by a tree process.

Reconstruction of the maximum-likelihood tree structure for single-cell transcriptomics datasets

As we focus on the analysis of scRNA-seq data, we will in the following presentation assume that the objects of our dataset are cells and their positions in high-dimensional space are their gene expression states, even though our methods apply much more generally.

For a given dataset D, consisting of measurements of the gene expression states of a collection of cells, we need to define a model that assigns a likelihood P(D∣T, t) to any tree, as defined by its topology T and branch lengths t. This likelihood should reflect the probability for the gene expression states of the cells at the leaves to have diverged along the branches of this tree. The optimal representation of the dataset then corresponds to the tree that maximizes this likelihood. The likelihood naturally decomposes into factors corresponding to probabilities for the movements along the branches of the tree and factors for the probabilities of the measurements given the positions of the cells at the leaves:

where x_i is the true position in gene expression space of node i and π(j) is the ancestor of node j in the topology T. Note that, in this definition, we are marginalizing over the unknown positions x_j of each node j. In the Supplementary Information (section B1), we describe the derivation of this expression in detail.

The first product in the likelihood (1) runs over leaf nodes and P(cell_i∣x_i) is the probability of the observed data for cell i assuming true gene expression state x_i. Bonsai assumes that the likelihoods P(cell_i∣x_i) can be approximated as a product of independent Gaussians for each component of the feature vector x_i, with independent standard deviation for each feature (which may also vary across cells). For example, for scRNA-seq data, we need to provide a most likely position in the high-dimensional gene expression space, with independent uncertainty estimates for each gene in each cell (Fig. 1a). In our applications, these are obtained directly from the normalization method Sanity²⁴. The second product runs over all branches in the tree and corresponds to the probabilities for the movement in gene expression space along each branch. In particular, if node j is connected by a branch of length t_j to ancestral node π(j), then P(x_j∣x_π(j),t_j) corresponds to the probability of moving from the gene expression state x_π(j) to state x_j of node j. Note that the positions of ancestral nodes can be interpreted as likely states of ancestor cells in developmental settings or as possible intermediate cell states in a cellular reprogramming context.

As we mentioned above, we know relatively little about the dynamical process that governs gene regulation and a main reason for gathering scRNA-seq data is to learn about this process in the first place. We take a Bayesian approach, assuming as little as possible about which movements in gene expression space are more or less likely and letting the data guide us to the most likely tree structure. This tree structure then provides information on the gene expression dynamics along its branches. More specifically, we make the following assumptions to determine the probabilities P(x_j∣x_π(j),t_j). First, we assume that movements in gene expression space can be described by a continuous Markov process. That is, we expect no discontinuous jumps in gene expression and we expect changes in state to only depend on the current gene expression state. Note that this still allows for a change in gene expression to be caused by a signal outside of the cell as long as the probability for this signal to occur can be predicted by the current gene expression state.

Second, although we recognize that the dynamics likely proceed along some manifold embedded in the high-dimensional gene expression space, we assume that this manifold is sufficiently smooth and sampled sufficiently densely that, at least for the relatively short distance along each branch, it is sufficiently flat that the Euclidean distance between the branch’s start and endpoint accurately approximates its length along the manifold. Lastly, we do not presume any preferred directionality for the gene expression dynamics; that is, a priori, movement in any direction is equally likely. In the Supplementary Information (section B1.4), we show, using classical results in Markov process theory^27,28, that this constrains our prior to being a homogeneous diffusion process so that the conditional probabilities P(x_j∣x_π(j), t_j) are given by Gaussians. This allows us to perform all the integrals over the positions x_i analytically, obtaining an analytical expression for the likelihood of any tree (Fig. 1b and Supplementary Information, section B2).

Very roughly, as discussed in the Methods, the log likelihood of a tree topology under our model is approximately proportional to the sum of the logarithms of its branch lengths. This can be intuitively understood as follows. Because each branch is a priori equally likely to point in any direction, the probability of the observed direction is inversely proportional to the surface area of a sphere in gene expression space with radius equal to the length of the branch. The likelihood of all branch directions is, thus, inversely proportional to the product of the surface areas of all corresponding spheres and its logarithm is, thus, proportional to the sum of the logarithms of the branch lengths.

Bonsai searches the space of all trees for the maximum-likelihood tree (Fig. 1c). As the number of possible trees grows superexponentially with the number of cells²⁹, it is computationally infeasible to perform an exhaustive search. Thus, we developed a search algorithm inspired by phylogenetic methods^30,31,32,33, which we describe in detail in the Methods and Supplementary Information (section C). Briefly, we start from a star tree that consists of one root to which all cell nodes are connected. From that tree, we iteratively add internal nodes so as to maximally increase the likelihood at each step and stop when the likelihood can no longer be increased by adding an internal node. After this, we continue to search the tree space for a better solution using so-called subtree prune and regrafting (SPR-moves) and nearest-neighbor interchanges (NNI-moves). In addition, to facilitate tree reconstructions on large datasets with millions of cells, we provide a backbone-based version of Bonsai in which we first reconstruct a Bonsai tree that connects only a subset of the cells, after which the remaining cells are sequentially placed on this backbone (Methods).

After finding the optimal tree, Bonsai reports its topology T and branch lengths t in the standard Newick format. To visualize and interactively explore the resulting trees, we also provide an application, called Bonsai-scout, which allows different visual layouts¹⁸ and a variety of downstream analyses.

On simulated data, Bonsai representations accurately recover cell-to-cell distances on all scales

To test whether Bonsai can recover differentiation trajectories, we first created six synthetic datasets where, in the ground truth, the gene expression patterns of cells had diverged along the branches of a tree. To ensure that these data realistically mimic real scRNA-seq data, we matched statistics of real scRNA-seq datasets including the distributions of total unique molecular identifier (UMI) counts and the means and variances of the expression levels across genes and cells. We also added both biological and measurement sampling noise in line with current understanding of the noise properties of scRNA-seq measurements^24,34,35 (detailed description in the Supplementary Information, section E). In the simplest dataset, gene expression patterns diverged along a simple binary tree with all branches of equal length and each gene’s expression diffusing independently along each branch of the tree. In the subsequent five datasets, we added increasing complexities, including unbalanced trees, random branch lengths and mimicking observed correlations in expression across genes so that the diffusion along the tree’s branches effectively occurs in a lower-dimensional subspace.

Figure 2a compares the visualizations of the most complex dataset, in which the ground truth is an unbalanced tree with varying branch lengths and in which the true expression levels are correlated across genes in a manner mimicking real data. Whereas Bonsai almost perfectly recovers the structure and differentiation trajectories in this dataset, the principal component analysis (PCA) and UMAP visualizations fail almost entirely on this task. Although UMAP clusters some groups of similar cells together, it completely fails to represent the relationships between or within these groups and tuning its parameters does not solve this problem. The PCA visualization only manages to capture the rough global structure along the two dimensions with the most variation, that is, separating the light-green clade of cells from the others along the first axis and separating the purple and red from the blue and orange clades along the second axis. Similar results were obtained for all the other simulated datasets (Supplementary Figs. 2–7, including results for the more recent tools PHATE³⁶ and DTNE³⁷). On these datasets, DTNE’s visualizations appear to capture a bit more of the large-scale structure in the data than either UMAP or PHATE (which capture virtually none). However, no other method comes anywhere near capturing the structural detail that Bonsai correctly reproduces.

Next, we asked to what extent the different visualizations accurately represent the distances between cells. As discussed above, for these datasets, we expect the distances between cells along the ground-truth tree to match the Euclidean distances between their gene expression states; Supplementary Fig. 8 confirms that this is indeed the case for all six datasets. Thus, we asked to what extent the cell-to-cell distances in the visualizations match the Euclidean distances between their ground-truth states in the high-dimensional gene expression space. In particular, we calculated, for each cell, the correlation between the true distances to all other cells and the corresponding distances in the representation. Figure 2b shows histograms of these correlations for the dataset of Fig. 2a, and Fig. 2c shows box plots of the distributions of correlations for Sanity, Bonsai, PCA and UMAP across the datasets (with Supplementary Fig. 9 showing box plots for PHATE and DTNE as well). We find that only Bonsai manages to accurately capture almost all pairwise distances across these datasets, which is remarkable given that these synthetic datasets are all heavily affected by the biological and measurement sampling noise, as in real scRNA-seq datasets. In contrast, all other visualization methods perform far more poorly, with Sanity typically performing second best, followed by PCA and DTNE. UMAP and especially PHATE perform very poorly in reproducing true pairwise distances. These results suggest that, in contrast to existing methods, Bonsai not only captures the trajectories along which cells have diverged in gene expression space but also accurately conserves the pairwise distances between the cells at all scales.

Lastly, if we use conventional preprocessing instead of preprocessing with Sanity²⁴ and explicitly incorporating error bars on the gene expression estimates, Bonsai’s accuracy is much reduced (Supplementary Fig. 10 and Methods).

The blessing of dimensionality

But what if the data do not derive from a tree structure? It is not clear whether Bonsai’s tree representations can also conserve cell-to-cell distances for datasets that were not generated along a tree structure. To test this, we created a non-treelike dataset by simulating seven ‘pseudobulk’ clusters with sizes ranging from 34 to 320 cells. Within each cluster, the cells were given almost identical ground-truth states but their raw data differ because the UMI counts of each cell had independent sampling noise. We find that the Bonsai tree accurately represents almost all cell-to-cell distances for this dataset as well (Fig. 2c and Supplementary Fig. 11). Thus, even though Bonsai’s tree optimizes the probability of the trajectories connecting the cells and these cells did not derive from a tree, the tree nonetheless accurately captures pairwise distances.

To test how generic this result is, we asked to what extent trees reconstructed on random scatters of points also conserve distances. In particular, we created datasets in which points were placed at random distances and in random directions from the origin in either 2, 10, 100, 1,000 or 10,000 dimensions. We ran Bonsai on these noiseless datasets and calculated the correlations between the true distances and those visualized on the tree (Supplementary Fig. 12).

We find that, while other visualization methods perform poorly and deteriorate quickly in performance as the dimensionality of the dataset increases, the accuracy of tree structures keeps improving and becomes essentially perfect from 1,000 dimensions onward. Thus, for random scatters of points, the higher the dimensionality, the better all pairwise distances are represented by the Bonsai tree. We next tested whether this observation extends to data with the structure of real scRNA-seq data. We created a ground-truth dataset by processing a dataset of cord-blood cells³⁸ using Sanity²⁴ and then treating the vectors of estimated log expression values as ground truth. From this noiseless ground-truth scatter in high dimensions, we created datasets of different dimensionality by projecting on the first 10, 50, 100, 500 or 1,000 PCA components. We then tested to what extent the cell-to-cell distances in the ground truth matched those of the Bonsai tree representations and other visualizations. We again find that, as the dimensionality of the dataset increases, the tree representations perform better and better in capturing true pairwise distances (Supplementary Fig. 13). Moreover, whereas the tree representation almost perfectly captures the cell-to-cell distances in the high-dimensional data, the distances in the representations of all other visualization methods hardly correlate at all with the ground truth.

Thus, we find that, even though Bonsai’s trees are optimized to find the most likely trajectories connecting the cells, these trees also conserve the pairwise distances between the cells, even if these cells did not derive from a tree. That is, rather than a curse, there is a blessing of dimensionality, whereby the higher the dimensionality of the data, the better all pairwise distances are represented by distances along the branches of the most likely tree connecting the cells.

Bonsai outperforms existing methods on nearest-neighbor identification

Figure 2b,c shows that, in addition to vastly outperforming other visualization methods, Bonsai also substantially outperforms the distance estimates of Sanity²⁴. As Sanity uses the same gene expression estimates and error bars and also explicitly incorporates the error bars in its estimates, this suggests that the tree reconstruction itself improves distance estimates by forcing the estimated distances to be consistent between nearby cells on the tree.

To test this directly, we started from the dataset with 100 randomly placed points in a 100-dimensional space of Supplementary Fig. 12 and then replaced each of the 100 points with a cloud of either 1, 2, 5, 10 or 20 noisy ‘measurements’ near the original point (creating datasets with 100, 200, 500, 1,000 and 2,000 points total). For each of these datasets, we then asked each method to estimate the true distances between the points. We find that Bonsai not only dramatically outperforms all other methods but also, in contrast to conventional embedding methods, Bonsai’s distance estimates improve as the number of measurements increases (Supplementary Fig. 14). This confirms that, by representing the similarities between the points on a tree, the estimated position of each point improves because of the presence of nearby points. Bonsai, thus, automatically accomplishes the type of regularization that data diffusion methods aim at by estimating states of individual cells by smoothing over the states of predicted nearest neighbors^39,40.

We next investigated the consequences of this effect for the identification of nearest-neighbor cells. Estimating nearest-neighbor cells is important because many scRNA-seq analysis methods use so-called k-nearest-neighbor graphs, where cells are connected to their k nearest neighbors in gene expression space. We previously showed that, because of the large measurement noise in scRNA-seq data, it is typically very challenging to correctly identify the nearest neighbors of each cell²⁴ and Sanity’s distance estimates substantially outcompete other methods on this task.

Here we compared Bonsai’s performance on identifying nearest neighbors to the performance of Sanity, PCA, UMAP, PHATE and DTNE on all seven synthetic datasets and find that Bonsai substantially outperforms these methods (Supplementary Fig. 15).

Bonsai scales to large cell atlases

Bonsai implements various mathematical and computational speed-ups (Methods and Supplementary Information, sections C and D) ensuring that its computation time and memory usage stay within acceptable limits. For example, the majority of scRNA-seq datasets have fewer than 30,000 cells (Suppl. Fig. 16) and take less than 1 day to process with the standard Bonsai algorithm (Fig. 3a, red dots). To allow scaling to very large datasets, we developed a backbone-based mode inspired by large phylogenetic reconstructions⁴¹ (Methods) that allow tree reconstructions for datasets of over 1 million cells (Fig. 3a, blue dots). In this mode, we first use the standard algorithm to reconstruct a tree on a subset of cells to serve as a backbone, after which the remaining cells are iteratively placed on this backbone. In Supplementary Fig. 17, we show that the results of this faster algorithm are consistent with the normal Bonsai algorithm; in Supplementary Fig. 18, we show that its computation time scales as C^1.15 for large datasets, while memory usage also remains within bounds (Supplementary Fig. 19).

As a demonstration of Bonsai’s computational efficiency, and as resources for the community, we reconstructed trees for a retina atlas with 1.2 million cells⁴² (Fig. 3b–d) and the well-known cell atlases of 70,000 mouse cells from the Tabula Muris dataset⁴ (Supplementary Fig. 20) and 100,000 epithelial cells from the Tabula Sapiens dataset⁵ (Supplementary Fig. 21).

Interactive visualization and exploratory downstream analysis with Bonsai-scout

In conventional two-dimensional embeddings of large datasets, points are inevitably packed together so closely that any visual analysis is influenced not only by the parameters of the visualization method but even by the order in which points are plotted. Fortunately, the hierarchical structure of trees provides a natural way to navigate and zoom in on different subtrees. Although there are tools for DNA sequence analysis that can visualize large trees⁴³, these are not suitable for exploring single-cell data. Therefore, we built Bonsai-scout, a customized application that facilitates exploring Bonsai trees, as we showcase in Fig. 3b–d and Supplementary Figs. 20 and 21 for the retina, mouse and human atlases.

Bonsai-scout supports different layouts of the tree, including the equal-daylight layout¹⁸ shown in Fig. 3 and Supplementary Fig. 20 and the more conventional dendrogram (Supplementary Fig. 21). When using the equal-daylight layout, the tree is shown on a hyperbolic disk with tunable curvature, where distances are naturally compressed toward the edges of the disk⁴⁴. The user can interactively control the visualization by recentering it so as to focus on one area of the tree (Fig. 3c), actively zooming in on a particular part, flipping branches, overlaying cell type annotations (Fig. 3b,c), showing the expression of any gene on the tree (Fig. 3d), etc.

Bonsai-scout also implements several methods for more sophisticated downstream analysis. One of the most commonly applied scRNA-seq analyses is to identify ‘cell type’ clusters of cells with clearly distinct expression patterns. As described in the Methods, Bonsai-scout implements both an unsupervised and a supervised method for identifying such clusters of cells. In the unsupervised method, the tree is divided into subtrees by iteratively cutting branches so as to minimize the sum of pairwise distances between cells in the subtrees. In the supervised clustering, Bonsai-scout takes a user-defined annotation of cells into clusters as input and divides the tree into subtrees so as to maximize the normalized mutual information (NMI) between the provided annotation and the subtrees. To test the performance of these clustering methods and compare to the performance of the popular Leiden and Louvain clustering methods, we collected all scRNA-seq datasets with fewer than 100,000 cells from the Bgee database⁴⁵ and calculated the match of the reference annotations with the clusterings obtained with each method (Supplementary Fig. 22). Bonsai-scout’s clustering results generally show a similar match with the reference annotations as the results of Leiden and Louvain clustering. This confirms that annotated cell types generally correspond to clades of the Bonsai tree.

Once clades of interest have been identified, Bonsai-scout can also be used to identify marker genes that distinguish between different clades of interest (Fig. 3d and Methods).

Lastly, one may ask what happens if the data lack any structure at all. Notoriously, clustering methods tend to return clusters even if there is no clear clustering structure in the data. To test Bonsai’s results on structureless data, we created random data by sampling a random integer between 0 and 1,000 for each gene and cell. Because of the high dimensionality of gene expression space, this creates data in which cells are almost perfectly equidistant from each other and a faithful representation should, thus, reflect this. In Supplementary Fig. 23, we show that Bonsai perfectly captures the equidistance in this dataset and reflects this absence of structure by approximately connecting all cells with long terminal branches to a single root. In contrast, other visualization methods infer a large variation in pairwise distances, placing some cells much closer to each other than others.

We make Bonsai-scout available both as a standalone program that users can run themselves and through a webserver where the Bonsai representations of several datasets are already available for exploration. Lastly, to facilitate Bonsai analysis, we have implemented an automated analysis pipeline online (https://bonsai.unibas.ch), where users can upload a UMI count table of any dataset of interest. The pipeline starts by running Cellstates²⁵ to identify groups of cells with statistically indistinguishable gene expression states, followed by Sanity²⁴ to obtain normalized gene expression levels (that is, log transcription quotients, LTQs) and their error bars. Then, these results are used to reconstruct the Bonsai tree, after which the user gets both flat file results and a link to a web interface where their results can be explored with Bonsai-scout.

Bonsai automatically identifies cells in precursor states or different developmental stages

What if an scRNA-seq dataset contains cells from different times or stages of a differentiation process? Because, per definition, Bonsai places every cell at a leaf of the tree and cells in precursor states rather correspond to internal nodes of the true differentiation tree, it is unclear whether Bonsai’s trees can readily identify cells at different temporal stages of differentiation. To test this, we took one of the simulated datasets where the cells at the leaves diverged along a binary tree but now also added a cell for each internal node of the tree to the dataset. Supplementary Fig. 24 shows that Bonsai connects each precursor cell by a short terminal branch to its corresponding internal node of the tree, so that the distance from the root correlates almost perfectly with each cell’s differentiation stage.

Next, we asked whether Bonsai can also successfully represent different development stages in a real developmental dataset. To test this, we ran Bonsai on a scRNA-seq dataset from Caenorhabditis elegans that combines data from different worms at different developmental stages and for which the authors provided an estimated developmental age of each cell⁴⁶. As shown in Supplementary Fig. 25, Bonsai detects multiple clades in the tree reflecting the start of different cell type lineages and the time ordering is conserved along the branches of each of these lineages. In particular, the distance of cells from the root in Bonsai’s tree correlates well with their annotated development stage.

Bonsai infers differentiation trajectories of blood cell types

For a more in-depth test of Bonsai’s ability to recover differentiation trajectories, we decided to focus on a dataset of 7,509 cord-blood mononuclear cells³⁸ because, arguably, the differentiation of blood cell types is one of the most studied mammalian differentiation processes. That is, there is considerable prior knowledge about the cell type differentiation hierarchy and on which genes can be used as cell type markers^47,48,49. In addition, this so-called CITE-seq dataset combines the measurement of gene expression (scRNA-seq) with measurements of a selection 13 cell surface markers that are often used for cell type classification, allowing independent validation of the cell types of individual cells and, thus, of the structures Bonsai infers on this scRNA-seq data.

In Fig. 4a, Bonsai’s representation is visualized as a dendrogram with cells colored on the basis of the cell type annotations from the original publication³⁸. We see that the Bonsai tree clearly shows a major split between a large myeloid clade comprising monocytes and dendritic cells and a clade with lymphoid cell types. The identity of the myeloid clade is independently confirmed by the abundance of the myeloid surface marker protein CD11c (Supplementary Fig. 26). Among the lymphoid cells, Bonsai also identifies distinct clades with known cell types; for example, we see a distinct clade of T cells, as confirmed by the abundance of the CD3 marker protein complex (Supplementary Fig. 26). We also observe that red blood cells and megakaryocytes each form distinct clades that are distant from all other cell types, consistent with their unique transcriptional profiles, dominated by hemoglobin expression in red blood cells and by platelet production in megakaryocytes.

To investigate in more detail to what extent Bonsai recovers known differentiation trajectories of the major blood cell types, we reduced the complexity of the data by first running Cellstates²⁵, which groups cells with statistically indistinguishable gene expression profiles. Cellstates not only reduces the number of states but also, by aggregating counts from cells in these ‘cell states’, we reduce the measurement noise, so that the preprocessing performed with Sanity and Bonsai becomes more accurate. Running Bonsai on the output of Cellstates results in a lineage tree of the major cell types (Fig. 4b) that is in remarkable agreement with the current view in the literature (Fig. 4c, adapted from a previous study⁵⁰). Indeed, starting from the hematopoietic stem cells at the root, erythrocytes branch off in the first branching, followed by the major split between myeloid and lymphoid cells. The lymphoid cells further separate into T, B and natural killer (NK) cells, while the myeloid cells split into dendritic cells and monocytes. The only major discrepancy between the Bonsai tree and accepted lineage relationships is the placement of the megakaryocytes, which Bonsai places in the myeloid branch as opposed to the branch with the erythrocytes.

It is instructive to compare the detailed structure and differentiation hierarchy that Bonsai recovers to the visualizations of the same data by other methods (Supplementary Fig. 27). UMAP, PHATE and DTNE group the cells into various connected blobs that reflect cell type annotation but it is unclear what the shapes and locations of these blobs tell us about the heterogeneity within and between cell types. Moreover, there is virtually no agreement in these features across the methods. For example, the DTNE visualization suggests that the hematopoietic stem cells (HSCs) are so heterogeneous that some are more distant from each other than their distance to any other lymphoid cell type (that is, B, T and NK cells). In contrast, the HSCs form a tight blob in the UMAP visualization. Furthermore, UMAP suggests HSCs are closest to a subset of erythrocytes and B cells, while PHATE suggests that they are far from most other cell types. Similarly, while Bonsai correctly infers that expression profiles of erythrocytes and megakaryocytes are particularly far from those of the other cell types, this is not clear at all in any of the other visualizations. Lastly, while Bonsai makes explicit predictions about the differentiation hierarchy relating the cell types that closely match the literature, the other visualization methods leave the differentiation hierarchy fully up to the imagination of the user.

Returning to Bonsai’s predictions, we note that, instead of a single clade of NK cells, the trees of Fig. 4a,b place NK cells in two groups, with one located in the lymphoid branch of the tree and a smaller group within the myeloid branch (Fig. 4d). While it is generally assumed that NK cells solely derive from lymphoid progenitors^{51,52,53,54,55}, in vitro studies have shown that NK cells can be generated from myeloid progenitors^56,57 and it has been debated whether such myeloid-derived NK cells also occur in the in vivo developmental trajectory of NK cells⁵⁸. Thus, Bonsai’s results strongly suggest that some NK cells indeed do derive from the myeloid branch.

As a first validation of this prediction, we used the measured abundances of surface marker proteins to confirm that, in addition to the myeloid marker (CD11c), these cells express the canonical NK cell surface markers CD56 and CD16, as well as NK marker genes at the mRNA level such as NKG7 (Fig. 4e). That is, these cells indeed carry known markers of both myeloid origin and NK cell function. To exclude the possibility that these cells are a technical artifact deriving from ‘doublets’ where two cells are inadvertently measured together, we reran our analysis after first removing doublets using Scrublet⁵⁹. We find that the tree structure is virtually unchanged and still contains the two NK groups (Supplementary Fig. 28). We further validated our finding by investigating whether myeloid-derived NK cells can also be found in other datasets. We analyzed a CITE-seq dataset containing about 30,000 human bone marrow cells⁶⁰ and found that myeloid-derived NK cells are indeed also present in this bone marrow dataset, albeit at substantially lower frequency, suggesting that myeloid-derived NK cells may be more abundant in cord blood (Supplementary Fig. 29). Together, these independent validations add substantial support to Bonsai’s prediction of a subtype of NK cells originating from myeloid progenitors in vivo.

Bonsai’s marker gene identification method allowed us to identify a substantial number of marker genes (Methods) whose expression clearly separates the lymphoid NK from myeloid NK cells (Fig. 4f and Supplementary Figs. 30 and 31). Looking at these marker genes, we may even speculate about potential functional differences between myeloid NK and lymphoid NK cells. It appears that the genes that are upregulated in the myeloid NK cells are more related to the recognition of ‘foreign’ molecular patterns, metal ion sequestration and release of granules with antimicrobial activity, while the genes that are up in the lymphoid NK cells are more involved in cytolysis (cell killing) and the destruction of virus-infected and tumor cells. Notably, markers of the lymphoid NK cells include major histocompatibility complex (MHC) class I genes, whereas markers of the myeloid NK cells include MHC class II genes. This is consistent with this proposed functional subdivision, as MHC class I molecules present intracellular antigens to CD8⁺ cytotoxic T cells, which can kill infected or cancerous cells, whereas MHC class II molecules primarily present extracellular antigens to CD4⁺ helper T cells. While these are just speculations at this point, they underscore that Bonsai not only predicts previously undescribed subtypes of cells but also identifies marker genes distinguishing these subtypes, suggesting specific biological hypotheses regarding the roles of these subtypes.

As we do not a priori know what complex high-dimensional structures may be hidden within single-cell omics datasets, there is an urgent need for exploratory analysis and visualization tools that accurately represent these structures without distortion. Unfortunately, for lack of a better alternative, most researchers continue to use visualization tools^9,10 that are well known to severely distort the true structure in the data and hallucinate structure that does not exist^11,12,13.

Motivated by a long history of representing high-dimensional objects using hierarchical structures and the knowledge that the gene expression states of cells in an organism have in fact diverged along the branches of a cell division tree, we here proposed to represent single-cell omics data using trees with cells at the leaves. Starting from a minimal number of assumptions, we derived a Bayesian model without any tunable parameters that finds the tree that maximizes the probability of the gene expression changes along the branches of the tree while taking into account the heterogeneous measurement noise properties of the data at the leaves.

We showed that the resulting Bonsai method successfully visualizes the structures in scRNA-seq datasets without distortion. On realistic simulated data, Bonsai recovers ground-truth differentiation trajectories with high accuracy, while other methods are unable to capture relationships beyond clustering similar cells together (Fig. 2 and Supplementary Figs. 2–11). Although Bonsai was designed to recover differentiation trajectories, we found that, because of an effect we dubbed the blessing of dimensionality, its tree representations also accurately conserve true distances in the high-dimensional gene expression space without distortion (Supplementary Figs. 12 and 13). Moreover, in contrast to other methods, Bonsai’s trees faithfully positions cells in precursor states (Supplementary Fig. 24) or different developmental stages (Supplementary Fig. 25). On a dataset of cord-blood cells, Bonsai not only recovers known lineage relationships with high accuracy but also discovers a subtype of natural killer cells deriving from the myeloid lineage that, as far as we are aware, has not previously been described in vivo (Fig. 4).

Notably, the similarly poor performance of other methods in capturing more global relationships across our real and simulated datasets is explained by the fact that these methods are all based on similar principles, that is, mainly aiming to conserve nearest-neighbor relationships²⁶. This also applies to existing methods for trajectory inference that are based on nearest-neighbor graphs or low-dimensional embeddings of the data^61,62,63, which are heavily affected by noise and have proven to be shaky foundations for further analysis^{11,12,13,64,65}. In fact, we showed that Bonsai also strongly outperforms other methods on the basic task of identifying nearest neighbors in the first place (Supplementary Fig. 15).

Bonsai’s tree representations of scRNA-seq data now enable many different downstream analyses that would otherwise be very challenging to implement. For example, by inferring the gene expression states at all internal nodes of the tree, Bonsai directly quantifies the expression changes that occurred along each branch of the tree. Indeed, it is well established in the field of evolutionary theory that, to understand evolutionary dynamics, one should focus on the changes that occur along the branches of the phylogeny⁶⁶. Similarly, while coregulation of genes is typically inferred by treating all cells as independent samples, it is more informative to study the correlation structure of gene expression changes along the branches of the tree because these reflect the differentiation trajectories of the cells and, thus, more directly reflect the actions of the gene regulatory circuitry. Bonsai’s tree representations can also be used to locate major differentiation events or to systematically estimate densities of cells in gene expression space.

There are many avenues for future extensions of the Bonsai analysis. First, Bonsai assumes that Euclidean distance is an accurate metric for the relatively short distances between neighboring nodes on the tree and that, a priori, changes in any direction in the high-dimensional gene expression space are equally likely. However, because gene expression changes are to a large extent driven by changes in the activities of regulators, we not only expect gene expression states to be constrained to a lower-dimensional manifold but also expect gene expression changes to exhibit a nontrivial correlation structure. Indeed, on real data, Bonsai readily infers substantial expression correlations between genes (Supplementary Fig. 37). In addition, although Bonsai still performs robustly on simulated data that include realistic gene expression correlations, we found that its accuracy decreases when such correlations exist (Fig. 2c). Bonsai’s accuracy could, thus, likely be improved by learning such a correlation structure from the data to incorporate it into its probability model for movement in gene expression space.

Second, many single-cell omics studies involve batches of cells that derive from different biological samples (for example, different animals) and there may be systematic differences in the gene expression states of cells from different batches that complicate the interpretation of the results. However, the treatment of such ‘batch effects’ in the analysis will be highly dependent on the specific experimental design and motivating scientific questions. For example, in one study, systematic differences across samples from different individuals might be a nuisance that needs to be corrected for, whereas, in another study, the differences across individuals might be the main signal of interest. Currently Bonsai implements a simple batch correction that corrects for systematic differences in average gene expression levels across batches (Methods) but more complex batch correction procedures are a possible avenue for extension of Bonsai.

Third, Bonsai could be applied to multimodal datasets. For example, it is currently possible to obtain both transcriptomic (scRNA-seq) and epigenomic (scATAC-seq) readouts from the same cells⁶⁷; however, because scATAC-seq data are even more sparse than scRNA-seq, such data are typically analyzed by first grouping cells into so-called pseudobulk samples. Bonsai provides a principled way of doing this by using the clades of the scRNA-seq tree to define pseudobulk samples. Moreover, with suitable preprocessing, Bonsai could be extended to build trees directly from the scATAC-seq data, allowing direct comparison of Bonsai trees built from transcriptomic and epigenomic cell states.

Fourth, it is increasingly appreciated that having lineage information at the single-cell level is crucial for understanding gene expression dynamics^68,69 and there are now various experimental techniques that measure single-cell transcriptomes while simultaneously recording lineage histories by creating random mutations on genetically encoded barcodes^{70,71,72,73,74,75}. An exciting prospect is to extend Bonsai to combine gene expression information with lineage-tracing data by extending the likelihood function to include the observed mutations in the lineage-tracing barcodes. This would not only allow more accurate lineage reconstructions but also allow systematic study of the extent to which gene expression and cell division histories align.

Last, in recent years, tools such as UMAP have become popular for visualizing not only scRNA-seq and scATAC-seq data but also many other types of biological data, including flow cytometry¹³, metagenomics⁷⁶, population genetics⁷⁷, neural spiking data⁷⁸ and even image and text categorization⁷⁹. Although our results focused on scRNA-seq datasets, the blessing of dimensionality that we identified above strongly suggests that, in principle, Bonsai’s trees should allow accurate visualization of any high-dimensional dataset with continuous features. Moreover, like the scRNA-seq data that we focused on here, many of these data types suffer from heterogeneous measurement noise properties that Bonsai is especially equipped to deal with.

To illustrate that Bonsai can robustly identify the structure of any high-dimensional dataset, we applied it to a dataset with statistics of professional football players (https://www.kaggle.com/datasets/vivovinco/20212022-football-player-stats) and found that the structure of Bonsai’s tree matches our intuitive understanding of the roles of different players (Supplementary Information, section H, and Supplementary Figs. 32–36). Thus, we expect that, for a wide variety of high-dimensional data, Bonsai’s tree representations will be able to much more accurately capture their structure, thereby providing a powerful method of exploratory data analysis across a wide variety of scientific fields.

In this section, we give an overview of the Bonsai method, its implementation and how specific datasets were analyzed. More detailed mathematical derivations underlying Bonsai and details of its computational implementation are provided in the Supplementary Information.

Recursive calculation of the tree likelihood

The foundation of Bonsai is a likelihood model for the observed data of each cell i given a tree topology T and a set of branch lengths t. This likelihood takes the form

where the integrals are over the unknown true positions x_j of each node j in gene expression space. As described in the Supplementary Information (section B), a crucial ingredient to the tractability of this model is that all factors in the likelihood have a Gaussian form that also factorizes over the genes. First, the likelihood of the observed data cell_i given the associated node’s true position x_i is given by

where μ_gi is the estimated cell position from the data and σ_gi is the standard deviation on this estimate for each gene g. Note that the size of the noise σ_gi may generally vary both across genes and cells. Second, the likelihood of the gene expression change along the tree branch leading from node π(i) to i is given by

where x_gi and x_gπ(i) give the coordinates of the two nodes, v_g is the estimated variance in the expression of gene g and t_i is the length of the branch. Note that the expected square of the change in gene expression of each gene g is proportional to the branch length t_i and the overall variance v_g of the gene.

As described in the Supplementary Information (section B2), we found that, for any tree topology, this likelihood can be efficiently calculated recursively, giving a continuous version of Felsenstein’s pruning algorithm²². For each node, say j, we can express the likelihood of all the data downstream of this node in the tree as a Gaussian function with effective positions \({\overline{x}}_{gj}\) and associated standard deviation \({\overline{\sigma }}_{gj}\) for each gene g. We can, thus, replace the likelihood contribution of the entire subtree downstream of node j by the likelihood of an effective leaf. This means that we can recursively simplify the tree by pruning off subtrees and summarizing their likelihoods as effective leaf contributions. Eventually, we get the likelihood of all data as a function of the root’s position P(D∣T, x_root), where we have marginalized over all other node positions. Importantly, because we assume a uniform prior on the root’s position, this expression is also directly proportional to the posterior distribution on the root’s position P(x_root∣D, T), which we exploit below. However, to get the likelihood of the tree T, we additionally marginalize over the root’s position.

Importantly, thanks to this pruning procedure, all integrals in Eq. (2) can be performed analytically, allowing the likelihood P(D∣T, t) to be calculated efficiently.

Intuitive approximation of the likelihood function

It is helpful to develop an intuition regarding how the likelihood of a tree topology depends on the distances along the branches of the tree. To do this, we calculate a log likelihood of a tree topology T by considering the simplified situation where there is no measurement noise, that is, the position x_i of each cell i at the leaves is known precisely. In addition, instead of marginalizing over all the positions of the internal nodes and optimizing the branch lengths, we optimize both with respect to all internal node positions and all branch lengths.

We then find that the log likelihood of a tree topology T is simply given (up to an additive constant) by a sum over the logarithms of the squared distances along the branches of the tree, that is,

where G is the number of genes, π(j) is the ancestor node of node j and the positions x_j at the internal nodes are set so as to maximize this likelihood. This occurs when the position x_j of each internal node j is equal to the weighted average of the nodes that it is connected to, with weights given by the inverse squared distances to each of the nodes. That is, the optimal internal node positions x_j obey

where the squared distances along the branches are represented by

and C(j) denotes the set of nodes that are connected to node j.

Although this simple expression for the likelihood of a topology only holds when there is no measurement noise on the leaves and when we optimize not only the branch lengths but also the positions of the internal nodes (rather than marginalizing over them), it illustrates that, under our model, the log likelihood does not correspond to a sum of distances or squared distances along the branches but is rather approximately equal to the sum of the logarithms of the distances along the branches.

Why the sum of the logarithms? This can be understood as follows. Optimizing the branch lengths and positions of the internal nodes causes the length of each branch to perfectly match the squared distance between its child and parent node. The only part of the data that remains unexplained by the tree is the direction of the change along each branch. As all directions are a priori equally likely, the probability of the observed direction is the inverse of the surface area of a sphere in G-dimensional space with radius equal to the length of the branch, which is proportional to the branch length to the power G-1. The probability of the entire tree is, thus, the inverse of the product of the surface areas associated with each branch and its logarithm is proportional to minus the sum of the logarithms of the branch lengths.

Obtaining posteriors over the positions of the internal nodes

In the Supplementary Information (section B1.5), we show that the tree likelihood is independent of which node we choose as the root of the tree. Therefore, by picking a node of interest i as the root and marginalizing over all other node positions as outlined above, we can get a posterior distribution P(x_i∣D, T, t) for the position x_i of each node i.

This posterior gives the probability of the internal node being in a certain gene expression state, conditioned on the data and the tree topology with corresponding branch lengths. These posteriors all have Gaussian form and Bonsai reports their mean and standard deviation, which can then be used for downstream analysis, such as in the Bonsai-scout tool for data visualization and marker gene detection.

An outline of Bonsai’s tree-search algorithm

We briefly outline how we search the large space of possible trees for the tree topology and branch lengths that jointly maximize the likelihood P(D∣T, t) (details in the Supplementary Information, section C). The main steps in Bonsai’s tree-search algorithm are as follows:

1. 1.

   Start from a star tree with optimized branch lengths.

2. 2.

   Iteratively add internal nodes to maximally increase the likelihood.

3. 3.

   Resolve polytomies.

4. 4.

   Reoptimize the branch lengths.

5. 5.

   Locally search for higher-likelihood trees by SPR-moves.

6. 6.

   Further increase the likelihood by interchanging nearest neighbors.

7. 7.

   Perform a final reoptimization of the branch lengths.

To start, we create a tree where each data point has an associated leaf node and these leaf nodes are all connected to a root node (that is, a star tree). For this star tree, we can efficiently optimize the branch lengths, after which we iteratively add internal nodes to maximally increase the likelihood at each step.

Importantly, it is also possible to incorporate prior information by starting with a different initial tree. For example, we recommend for scRNA-seq data to first calculate Cellstates²⁵ to group cells that are statistically indistinguishable. One can then start from a tree in which, for each cell state, the cells in this cell state are connected to an ancestor node and each of these ancestors is in turn connected to the root node, thereby initially forcing the cells within each cell state to form a clade in the tree.

From the initial tree, we iteratively pick a pair of leaves and add an internal node upstream of these two leaves so as to maximally increase the likelihood of the tree (when the branch lengths from the internal node to the two leaves and to the root are optimized). After adding the internal node, we use the pruning procedure described above, replacing the internal node with an effective leaf with a corresponding effective position and effective uncertainty (explanatory illustration in Supplementary Fig. 38). After this, the tree is again a star tree and we can iterate the procedure of adding an optimal internal node until we either have a fully resolved tree or can no longer add any internal node directly downstream of the root that increases the likelihood.

Polytomies can occur either because we started from an initial tree with polytomies or can be created when an optimal branch length is zero. For example, if an ancestor a is placed upstream of node i and leaf l, but with branch length zero from a to i, then this creates a polytomy with branches to the children of i and the leaf l all diverging from a. We attempt to increase the likelihood further by resolving such polytomies. We do this by effectively treating the internal node with the polytomy as a new root and then going over the above iterative procedure of adding internal nodes that increase the likelihood most.

Although these greedy procedures generally lead to trees with high likelihood, we found that, in most cases, local rearrangements can still further increase the likelihood. We search for such local rearrangements by performing SPR and NNI moves. In an SPR-move, we can pick any node of the tree (apart for the root) and disconnect it and its downstream subtree from the tree (subtree pruning). Then, we search over the remaining tree to find the best node to reattach this subtree (regrafting). In the Supplementary Information (section C5), we describe the different schemes that we use for selecting the candidate subtree to prune and for searching the best regrafting position. Every NNI-move starts by picking an edge, say from node a to b, and listing all the nearest-neighbor nodes, that is, nodes that are connected to either a or b, and then considering all ways of reconnecting these nodes to a and b. Bonsai first performs a random NNI-phase in which random edges are chosen and reconnections are sampled weighted by their corresponding change in likelihood, which is followed by a greedy NNI-phase on which the NNI-move that most increases the likelihood is performed deterministically.

Optional parameters

Bonsai does not require the user to specify any parameters; however, if desired, there are a few optional variants of its default behavior that the user can specify. First, for typical scRNA-seq datasets, there are many genes that are sampled so sparsely that the error bars on their expression measurements are larger than the true variation in their expression values and we found that including these genes in Bonsai’s analysis is more likely to decrease than increase the performance (likely because the Gaussian approximation of the measurement noise for these low signal genes is inaccurate).

From the raw UMI counts, Sanity calculates estimated gene expression values \({x}_{gi}^{* }\) for each gene g in each cell i and associated error bars ϵ_gi. These expression estimate \({x}_{gi}^{* }\) correspond to LTQs, that is, the expected logarithm of the fraction of the mRNA pool in cell i that corresponds to mRNAs of gene g. From these estimates, Bonsai calculates an average signal-to-noise ratio S_g for each gene as

where C is the number of cells and x_g is the estimated average LTQ of gene g across cells. For data other than scRNA-seq data processed by Sanity, Bonsai estimates the signal-to-noise ratio of each feature from the input data as described in the Supplementary Information (section B1.3). By default, only genes with S_g ≥ 1 are included in the analysis but the user can change this threshold value.

Second, by default, Bonsai uses a prior where the amount of diffusion in dimension g is assumed proportional to the variance v_g in the expression of gene g (Supplementary Information, section B1.4). If desired, the user can specify to not scale the prior by these variances v_g.

Correcting for differences in average gene expression across batches

Biological datasets often comprise data originating from different ‘batches’, for example, from different experiments, different animals or even different labs, and such batches can show various systematic differences. As explained in the Discussion, whether and how one wants to correct for these batch effects often depends on the experimental design and also on the biological questions that one is seeking to answer, whereby it is impossible to provide a general batch correction method that will be appropriate in all imaginable situations. However, we implemented an adapted version of scRNA-seq preprocessing that corrects for systematic differences in the average expression levels of genes across batches. This batch correction procedure is appropriate when fold changes across cells in each batch are more meaningful than systematic changes in the average expression levels of genes across batches.

Specifically, starting from a matrix with raw mRNA counts, we provide a script that splits the matrix into a count matrix for each batch and runs Sanity separately for each batch. From this, we obtain an estimated mean expression value for each gene in each batch, \({\mu }_{g}^{(b)}\), and log fold changes, δ_gc, away from this batch mean for each gene and cell. A second script then aggregates the data from different batches by concatenating the log fold changes, δ_gc, while the batch means \({\mu }_{g}^{(b)}\) are replaced by overall dataset means μ_g. Lastly, note that, to correctly aggregate these log fold changes, it is important that we first correct for the Sanity prior per batch because Bonsai expects means and standard deviations that describe the likelihood function rather than the posterior (as described for the normal Bonsai procedure in the Supplementary Information, section B1.2, and Supplementary Eq. (5)). As this prior correction is now already complete, Bonsai should be run with the argument --‘input_is_sanity_output False’. To ensure that this is applied correctly, our batch correction script already generates a configurations file that can be used for running Bonsai on the batch-corrected data. All necessary scripts can be found in the subdirectory ‘Bonsai-data-representation/optional_preprocessing/batch_correction’ of our repository (https://doi.org/10.5281/zenodo.20370956)²³.

Backbone-based Bonsai

To facilitate reconstructing Bonsai trees on the largest datasets, we implemented a backbone-based version of Bonsai that trades off a slight decrease in accuracy for much faster run times. The main steps of the backbone-based Bonsai algorithm are as follows:

1. 1.

   Preprocess the input for all cells.

2. 2.

   Reconstruct a backbone tree on a random subset of cells.

3. 3.

   Place the remaining cells on the backbone tree.

4. 4.

   Refine the tree with all cells.

In the first step, preprocessing of the full dataset is performed exactly as in the standard Bonsai algorithm. Next, we select a random subset of cells and use the standard Bonsai algorithm to create a backbone tree on the basis of this selected subset of cells. In the third step, we iteratively add all remaining cells to the backbone tree, using a tree placement algorithm that is also used for finding good regrafting positions in the SPR-moves, as described in detail in the Supplementary Information (section B4).

Roughly, we first find a number of central nodes spread over the backbone tree that will be used as start points for searching the best tree placement; in Bonsai, the number of such start points is equal to the logarithm of the number of nodes in the tree. For each start point, we calculate an ‘attachment likelihood’, which is the resulting tree likelihood if a cell would be attached with an optimized branch length to that node. Then, we repeat this calculation for the neighboring nodes of the start point and we continue the search in the direction of all neighbors with attachment likelihood over a certain threshold. In this way, we search all paths as long as the attachment likelihood either increases or only decreases slightly, thereby using the backbone tree as a decision tree.

After this search, we attach the new cell to the node that provides the highest attachment likelihood. In addition, we check whether the likelihood can be further increased by attaching the new node to one of the connecting edges, rather than to the selected node itself. Because the backbone tree can change appreciably when many cells are added, we also perform a global reoptimization of the branch lengths and recalculation of the search start points after a fixed percentage of backbone tree growth.

Lastly, we take the tree that was generated by adding the remaining cells and use it as an initial tree in the standard Bonsai algorithm. In practice, as there are very few nodes in this final tree with polytomies, this mostly boils down to a final pass of SPR-moves, NNI-moves and optimization of the branch lengths.

Our implementation of backbone-based Bonsai includes a few parameters that affect the computational requirements and that users can tweak depending on the amount of computational resources available. First, the user can determine the number of cells in the initial backbone. We recommend taking a large enough subset of cells such that the backbone is representative for the structure in the full tree (default: 10,000 cells).

In addition, the user can choose to not add all remaining cells in one step but, rather, only grow the backbone tree by a certain factor. For example, if there were 250,000 cells in total, the user can choose to start with a backbone of 10,000 cells, only add 40,000 cells and then proceed to the final refinement step. The refined 50,000-cell tree can then be used as a backbone for adding the remaining 200,000 cells. This two step procedure comes at a computational cost but might allow the algorithm to find a higher-likelihood solution.

Finding clusters in the tree

Bonsai-scout, the interactive app that can be used for exploratory analysis and visualization of Bonsai trees, offers both an unsupervised and supervised method for annotating clusters of cells on the basis of the clade structure of the tree. The unsupervised method is intended for users that do not already have their own annotations of their cells, while the supervised method identifies the clades in the tree that best match a user-provided annotation.

In the unsupervised method, branches are iteratively cut so as to greedily minimize the sum of pairwise distances between all leaves in the created subtrees. More specifically, for any tree, we can calculate the sum of pairwise distances between all leaves along the branches of the tree. Note that this sum differs from the sum of all branch lengths because branches in the middle of the tree are traversed by many more leaf-to-leaf paths than branches near the leaves. Starting from the original tree, we create two subtrees by cutting the branch that minimizes the combined sum of pairwise distances between leaves in the two resulting subtrees. Given these new trees, we can now cut another branch to again minimize the pairwise distances in the three resulting subtrees, and so on. To create n clusters, one repeats this cutting procedure n − 1 times.

In the supervised method, we search for a way of cutting the tree into subtrees such that the NMI between subtree membership and the labels of the user-provided annotation is maximized. Specifically, we search over all possible subtree sets by combining ‘cutting’ moves, in which a branch of the Bonsai tree is cut to create an additional subtree, and ‘gluing’ moves, in which such a cut is reverted. At each step, we first randomly pick whether we do a cutting or a gluing move, after which we sample the branch to cut or glue in proportion to the change in clade entropy it causes. With this, we preferentially cut branches that create relatively large subtrees, rather than repeatedly cutting off small subtrees or single leaves, as we found that this leads to a more efficient search for the optimal subtree set. Given this selected move, we calculate the change in NMI that it causes: ΔNMI. If this change is positive, we always perform the move; if it is negative, we accept the move with probability \(\exp (\ \Delta \,\text{NMI}\,/T)\), where T is a randomness parameter that we slowly decrease during the search. Then, when this search converges, we perform a final greedy phase in which the moves that most increase the NMI are performed until no move further improves the NMI.

Choosing a root node

The final tree that Bonsai reports (in Newick format) is rooted. By default, we position the root on the branch that would be the first to be cut in the unsupervised clustering procedure described above.

As the tree’s likelihood is independent of the choice of a root node, this choice is somewhat arbitrary. Notably, the choice of root affects none of the results but it of course influences the visualization of the tree. When using our Bonsai-scout visualization tool, the root determines the initial centerpoint for circular layouts and it determines the leftmost point in a dendrogram layout. For the circular layout users can double-click to focus the visualization on any arbitrary point in the tree; for the dendrogram layout, the root can be repositioned interactively.

Detection of marker genes

Bonsai can also be used to detect marker genes or features that best distinguish the cells or objects from two clades of the tree. We used this feature to identify marker genes that distinguish between the two groups of NK cells in the blood cell dataset (Supplementary Figs. 30 and 31) and to annotate the different football player clusters (Supplementary Figs. 33–36).

Let C₁ and C₂ denote two clades of cells on the tree. We define marker genes as genes that either maximize or minimize the probability that, when picking a random pair of cells, c₁ ∈ C₁ and c₂ ∈ C₂, the gene expression \({x}_{g{c}_{1}}\) is higher than \({x}_{g{c}_{2}}\). Thus, dropping the gene indices, we want to calculate the following for each gene:

To calculate this, we use the posteriors over the positions of the cell-associated leaf nodes that Bonsai provides. Specifically, for cell c₁, the posterior \(P({x}_{{c}_{1}}| D,T,{\bf{t}})\) is a Gaussian distribution with mean \({\overline{\mu }}_{{c}_{1}}\) and variance \({\overline{\sigma }}_{{c}_{1}}^{\,2}\); a similar posterior applies to cell c₂. To obtain the probability \(P({x}_{{c}_{1}} > {x}_{{c}_{2}})\), we note that the probability distribution of the difference, \({x}_{{c}_{1}}-{x}_{{c}_{2}}\) is also a Gaussian with mean \({\overline{\mu }}_{{c}_{1}}-{\overline{\mu }}_{{c}_{2}}\) and variance \({\overline{\sigma}}_{{c}_{1}}^{\,2}+{\overline{\sigma}}_{{c}_{1}}^{\,2}\). Therefore, we can express the desired probability in terms of the cumulative distribution of these Gaussians, which can be efficiently calculated using the following error function:

Marker genes are those for which M is close to 1 or 0. We provide a supplementary tool with the Bonsai code that takes a file that specifies two groups of cells and calculates the marker scores M for all genes. This supplementary tool can be found in the Bonsai code repository under ‘downstream_analyses/calc_marker_genes.py’.

As the marker score calculation involves a sum over all pairs of cells, for large clades, the calculation can become too slow for ‘on the fly’ exploratory analysis in Bonsai-scout. Therefore, in Bonsai-scout, the marker gene calculation is approximated for large clades by ignoring the uncertainties on the cell posteriors. In that case, the marker score reduces to

where Θ(•) is the Heaviside function, which is 1 for values larger than zero and 0 otherwise. Now, let \({N}_{{C}_{1}}(x)\) be the number of cells in C₁ that have an expression value lower than x, \({x}_{{c}_{1}}\) be the expression level of cell c₁ and rk(c₁) be the rank of cell c₁ in the list of all cells from C₁ and C₂ sorted in increasing order of expression. We then get

which we used in the second step to reorder the sum over c₁ such that we get a simple arithmetic sequence. This approximation of the marker score can be used in the interactive exploration of the data in Bonsai-scout, after which the expression from Eq. (10) can be used to confirm or falsify specific hypotheses. Notably, Bonsai-scout also allows for downloading files that can be directly used as input for the Python script that we provide to obtain the full marker scores.

Details of the processing of the different datasets

Processing of simulated datasets

For Bonsai, we first processed the raw UMI counts of the simulated datasets using Sanity²⁴ and then ran Bonsai with the default signal-to-noise threshold of 1.0. For the PCA, UMAP, PHATE and DTNE visualizations, we followed the most commonly used preprocessing steps. The raw counts were log-transformed with a pseudocount, \({x}_{gc}=\log (\frac{{n}_{gc}}{{N}_{c}}\overline{N}+1)\), where N_c is the total UMI count for cell c and \(\bar{N}\) is the average total UMI count of all cells. These log-transformed data were then projected on the first two PCA components for creating the PCA visualization and on the first N PCA components for UMAP, PHATE and DTNE, where N was optimized for each method by performing a coarse parameter scan on each simulated dataset. The results do not qualitatively depend on the number of PCA components used. We ran UMAP with the parameters n_neighbors = 15, min_dist = 0.1, n_components = 2 and metric = ‘euclidean’; PHATE and DTNE were run with their default parameters.

Lastly, to assess the effect of replacing the default Sanity preprocessing with simpler preprocessing, we also ran Bonsai on standard log-transformed counts with highly variable gene selection (using the same log transformation as above); we refer to this variant as Bonsai_log1p.

Processing of the cord-blood dataset

We downloaded the FASTQ files stored in the Gene Expression Omnibus (GEO) under accession number GSE100866, after which we pseudoaligned the reads with Kallisto using the parameter -x ‘0,0,16:0,16,25:1,0,0’. With bustools, we collapsed the UMIs and error-corrected the barcodes using a whitelist. Using the R package emptyDrops, we removed empty cells and only retained cells with a total UMI count of at least 500. From the remaining cells, we decided to keep only the cells that were also selected in the original publication such that cell type annotation was available. Lastly, for each gene, we counted the number of observed transcripts (that is, UMIs) in each cell and then ran Sanity²⁴ and Cellstates²⁵ on these raw counts.

Bonsai was then run on the log fold changes inferred by Sanity (stored in ‘delta_vmax.txt’) as gene expression estimates and using the corresponding standard deviations of the posteriors (stored in ‘d_delta_vmax.txt’) as the error bars. We started Bonsai from an initial tree informed by the Cellstates clusters; that is, we connected all cells from each Cellstates cluster to a single ancestor, which was in turn connected to the root.

Starting from the raw UMI count table, we ran PCA, UMAP, PHATE and DTNE visualizations using the standard preprocessing steps described above.

For visualizing the surface protein expression, we used the normalized values as described and provided in the original publication³⁸.

Processing of the Tabula Muris data

Raw FASTQ files were downloaded from the GEO under accession number GSE109774. All processing was conducted in the same way as described above for the cord-blood cell data.

Processing of the Tabula Sapiens data

For the Tabula Sapiens dataset, we downloaded the matrix of UMI counts from the CZI database (https://cellxgene.cziscience.com/collections/e5f58829-1a66-40b5-a624-9046778e74f5). We ran Cellstates²⁵ until the >104,000 initial cells were reduced to <5,000 clusters of cells with statistically indistinguishable gene expression states and used these to define an initial tree for Bonsai. The other processing steps were identical to those described in the previous sections.

Processing of the football statistics dataset

We obtained the dataset with football player statistics online (https://www.kaggle.com/datasets/vivovinco/20212022-football-player-stats). To preprocess this data, we first dropped players with unreliable statistics by (1) removing one player with NaN values; (2) keeping only the 50% players that played the most minutes to avoid outliers; and (3) removing 12 players that switched teams halfway through the season (unfortunately also dropping Cristiano Ronaldo). As the measured features are on vastly different scales, we also normalized them by first subtracting the mean for every feature and then rescaling such that the variance of each feature was equal to 1. After this, we projected the data on the first 50 PCA components. In this way, we reduced the effect of redundant statistics. The number of principal components used did not strongly affect the results. Then, because Bonsai expects uncertainty estimates on the input data but this information is not available for this dataset, we used an artificial small standard deviation of 10⁻³ for all features and players. The exact value chosen does not affect the Bonsai results either.

For the analysis of the Bonsai tree, we decided to use 12 clusters as they still gave interpretable groups, although we do not exclude the possibility that there is useful information in more fine-grained clustering. Outlier players were detected by first rooting the tree on the edge connecting the goal keepers to all other players. We then defined the 13 outliers by checking which players were furthest removed from the root. These outliers are highlighted in Supplementary Fig. 32.

An automated pipeline for Bonsai analysis of scRNA-seq data

We implemented an automated pipeline for scRNA-seq processing as a webserver (http://bonsai.unibas.ch). Here, users can upload a matrix with the raw mRNA counts per gene (rows) and cell (columns). In addition, users can upload annotations that they may already have for the cells and that can be visualized on the tree. All analysis is then performed automatically, which includes running Cellstates, Sanity and Bonsai. After all computations have finished, users get an email with a link to the Bonsai-scout visualization of their data and a link to a download page that contains flat files with all results.

Input and output of the Bonsai code

When using the Bonsai code independent of the automated pipeline, users should take note that Bonsai requires that the data are preprocessed such that the likelihood of the measurements is reasonably approximated by a multivariate Gaussian with means μ_gc, standard deviation σ_gc and negligible covariances (as discussed in the Supplementary Information, section B1). Bonsai requires the mean and standard deviation as simple feature-by-object matrices. For scRNA-seq data, the data can best be preprocessed using Sanity²⁴, after which Bonsai should just be pointed towards the full Sanity output folder.

The output that Bonsai provides contains the following:

• A Newick string with the tree topology and its branch lengths.

• Two files that describe the tree in a more human-readable format.

• Two Numpy binary files that describe the means and variances of the posteriors that Bonsai infers for the position of each node, that is, for all cells at the leaves and for the internal nodes.

• A .json file with metadata on the dataset, for example, containing the cell identifiers, gene identifiers and the inferred gene variances, as well as paths to where the original data were read from.

This output from Bonsai can be used as input for Bonsai-scout, for which we first run a preprocessing script that produces an .hdf file containing all necessary data and a .json file containing the initial settings for the tree visualization. As this .json file is human readable and editable, it is possible to change these initial settings by hand (for example, for customizing a color map). These two files are the only necessary files for running Bonsai-scout.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.