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

$$P(D| T,{\bf{t}})=\int\cdots \int\left(\prod _{i\in \,\text{cells}\,}P({{\rm{cell}}}_{i}\,| \,{{\bf{x}}}_{i})\times \prod _{{j\in\,\text{nodes}}\atop{j\ne \text{root}\,}}P({{\bf{x}}}_{j}\,|\,{{\bf{x}}}_{\pi (j)},{{t}}_{j})\right){\rm{d}}{{\bf{x}}}_{1}\cdots {\rm{d}}{{\bf{x}}}_{N},$$

(2)

where the integrals are over the unknown true positions xj 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 celli given the associated node’s true position xi is given by

$$P({\rm{cell}}_{i}| {{\bf{x}}}_{i})=\prod_{g}P({\rm{cell}}_{i}| {x}_{gi})=\prod_{g}\frac{1}{\sqrt{2\uppi}{\sigma }_{gi}}\exp \left(-\frac{{({x}_{gi}-{\mu}_{gi})}^{2}}{2{\sigma}_{gi}^{2}}\right),$$

(3)

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

$$P({{\bf{x}}}_{i}| {{\bf{x}}}_{\pi (i)},{t}_{i})=\prod _{g}\frac{1}{\sqrt{2\uppi {v}_{g}{t}_{i}}}\exp \left(-\frac{{({x}_{gi}-{x}_{g\pi (i)})}^{2}}{2{v}_{g}{t}_{i}}\right),$$

(4)

where xgi and xgπ(i) give the coordinates of the two nodes, vg is the estimated variance in the expression of gene g and ti 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 ti and the overall variance vg 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 algorithm22. 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, xroot), 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(xroot∣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 xi 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,

$$L(T)={\rm{cons.}}-\frac{G}{2}\sum _{j\ne {\rm{root}}}\log \left[\sum _{g}{\left({x}_{gj}-{x}_{g\pi (j)}\right)}^{2}\right],$$

(5)

where G is the number of genes, π(j) is the ancestor node of node j and the positions xj at the internal nodes are set so as to maximize this likelihood. This occurs when the position xj 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 xj obey

$${{\bf{x}}}_{j}=\frac{\sum _{k\in C(j)}\frac{{{\bf{x}}}_{k}}{{\Delta }_{jk}^{2}}}{\sum _{k\in C(j)}\frac{1}{{\Delta }_{jk}^{2}}},$$

(6)

where the squared distances along the branches are represented by

$${\Delta }_{jk}^{2}=\sum _{g}{\left({x}_{gj}-{x}_{gk}\right)}^{2},$$

(7)

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(xi∣D, T, t) for the position xi 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.

Start from a star tree with optimized branch lengths.

2.

Iteratively add internal nodes to maximally increase the likelihood.

3.

Resolve polytomies.

4.

Reoptimize the branch lengths.

5.

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

6.

Further increase the likelihood by interchanging nearest neighbors.

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 Cellstates25 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 Sg for each gene as

$${S}_{g}=\frac{1}{C}\sum _{i}\frac{{({x}_{gi}^{* }-{x}_{g})}^{2}}{{\epsilon }_{gi}^{2}},$$

(8)

where C is the number of cells and xg 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 Sg ≥ 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 vg 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 vg.

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)23.

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.

Preprocess the input for all cells.

2.

Reconstruct a backbone tree on a random subset of cells.

3.

Place the remaining cells on the backbone tree.

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 C1 and C2 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, c1 ∈ C1 and c2 ∈ C2, 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:

$$M=\frac{1}{| {C}_{1}| | {C}_{2}| }\sum _{{c}_{1}\in {C}_{1}}\sum _{{c}_{2}\in {C}_{2}}P({x}_{{c}_{1}} > {x}_{{c}_{2}}).$$

(9)

To calculate this, we use the posteriors over the positions of the cell-associated leaf nodes that Bonsai provides. Specifically, for cell c1, 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 c2. 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:

$$M=\frac{1}{| {C}_{1}||{C}_{2}|}\sum _{{{c}_{1}\in {C}_{1}}\atop{{c}_{2}\in {C}_{2}}}\left(\frac{1}{2}+\frac{1}{2}{\rm{erf}}\left(\frac{({\overline{\mu }}_{{c}_{1}}-{\overline{\mu }}_{{c}_{2}})}{\sqrt{2({\overline{\sigma }}_{{c}_{1}}^{\,2}+{\overline{\sigma }}_{{c}_{1}}^{\,2})}}\right)\right).$$

(10)

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

$$\tilde{M}=P\left({x}_{{c}_{1}} > {x}_{{c}_{2}}| {c}_{1}\in {C}_{1},{c}_{2}\in {C}_{2}\right)=\frac{1}{| {C}_{1}| | {C}_{2}| }\sum _{{c}_{1}\in {C}_{1}}\sum _{{c}_{2}\in {C}_{2}}{{\varTheta}} ({x}_{{c}_{1}}-{x}_{{c}_{2}}),$$

(11)

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 C1 that have an expression value lower than x, \({x}_{{c}_{1}}\) be the expression level of cell c1 and rk(c1) be the rank of cell c1 in the list of all cells from C1 and C2 sorted in increasing order of expression. We then get

$$\tilde{M}=\frac{1}{| {C}_{1}| | {C}_{2}| }\sum _{{c}_{1}\in {C}_{1}}\left(\,\text{rk}({c}_{1})-{N}_{{C}_{1}}({x}_{{c}_{1}})\right)=\frac{1}{| {C}_{1}| | {C}_{2}| }\left(\sum _{{c}_{1}\in {C}_{1}}\,\text{rk}({c}_{1})\right)-\frac{(| {C}_{1}| -1)}{2| {C}_{2}| },$$

(12)

which we used in the second step to reorder the sum over c1 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 datasetsProcessing of simulated datasets

For Bonsai, we first processed the raw UMI counts of the simulated datasets using Sanity24 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 Nc 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 Sanity24 and Cellstates25 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 publication38.

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 Cellstates25 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−3 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 Sanity24, 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.