Treefit - Getting started

This document is a user-friendly manual on how to use the Treefit package, which is the first toolkit for quantitative trajectory inference using single-cell gene expression data. In this tutorial, we demonstrate how to generate and analyze two kinds of toy data with the aim to help get familiar with the practical workflow of Treefit.

1. Generating toy data

While the Treefit package has been developed to help biologists who wish to perform trajectory inference from single-cell gene expression data, Treefit can also be viewed as a toolkit to generate and analyze a point cloud in d-dimensional Euclidean space (i.e., simulated gene expression data).

Treefit provides some useful functions to generate artificial datasets. For example, as we will now demonstrate, the function treefit.data.generate_2d_n_arms_star_data() creates data that approximately fit a star tree with the number of arms or branches; the term star means a tree that looks like the star symbol “*”; for example, the alphabet letters Y and X can be viewed as star trees that have three and four arms, respectively.

The rows and columns of the generated data correspond to data points (i.e., n single cells) and their features (i.e., expression values of d different genes), respectively. The Treefit package can be used to analyze both raw count data and normalized expression data, but regarding the production of toy data, it is meant to be used to generate continuous data like normalized expression values.

Importantly, we can generate data with a desired level of noise by changing the value of the fatness parameter of this function. For example, if you set the fatness parameter to 0.0 then you will get precisely tree-like data without noise. By contrast, setting the fatness to 1.0 gives very noisy data that are no longer tree-like. In this tutorial, we deal with two types of toy data whose fatness values are 0.1 and 0.8, respectively. We note that Treefit can be used to generate and analyze high dimensional datasets but we focus on generating 2-dimensional data to make things as simple as possible in this introductory tutorial.

1.1. Tree-like toy data

Let us first generate 2-dimensional tree-like data that contain 500 data points and reasonably fit a star tree with three arms. We can create such data and draw a scatter plot of them (Figure 1) simply by executing the following two lines of code:

In [1]: import matplotlib.pyplot as plt

In [2]: import treefit

In [3]: star_tree_like = treefit.data.generate_2d_n_arms_star_data(500, 3, 0.1)

In [4]: plt.figure();

In [5]: plt.scatter(star_tree_like[:, 0],
   ...:             star_tree_like[:, 1]);
   ...: 
_images/generate_tree_like_data.png

Figure 1. A scatter plot of tree-like toy data generated by Treefit

1.2. Less tree-like, noisy toy data

Similarly, we can generate noisy data that do not look so tree-like as the previous one. In this tutorial, we simply change the value of the fatness parameter from 0.1 to 0.8 in order to obtain non-tree-like data. The following two lines of code yields a scatter plot shown in Figure 2:

In [6]: star_less_tree_like = treefit.data.generate_2d_n_arms_star_data(500, 3, 0.8)

In [7]: plt.figure();

In [8]: plt.scatter(star_less_tree_like[:, 0],
   ...:             star_less_tree_like[:, 1]);
   ...: 
_images/generate_less_tree_like_data.png

Figure 2. A scatter plot of noisy toy data generated by Treefit

2. Analyzing the data

Having generated two datasets whose tree-likeness are different, we can analyze each of them. The Treefit package allows us to estimate how well the data can be explained by tree models and to predict how many “principal paths” there are in the best-fit tree. As shown in Figure 1, the points in the first data clearly form a shape like the letter “Y” and so the data are considered to fit a star tree with three arms very well, whereas Figure 2 indicates that the second data are no longer very tree-like because of the high level of added noise. Our goal in this tutorial is to reproduce these conclusions by using Treefit.

2.1. Tree-like toy data

Let us estimate the goodness-of-fit between tree models and the first toy data. This can be done by using treefit.treefit() as follows. The name parameter is optional but we should specify it, if possible. Because it’s useful to identify the estimation.

In [9]: fit_tree_like = treefit.treefit({'expression': star_tree_like},
   ...:                                 name='tree-like')
   ...: 

In [10]: print(fit_tree_like)
treefit.fit.Fit: tree-like
max_cca_distance:
     p      mean  standard_deviation
0    1  0.525533            0.256161
1    2  0.056173            0.012468
2    3  0.048170            0.006462
3    4  0.046294            0.006131
4    5  0.043460            0.005851
5    6  0.042430            0.005762
6    7  0.040809            0.005367
7    8  0.038870            0.005620
8    9  0.037796            0.005563
9   10  0.036458            0.005144
10  11  0.035436            0.004998
11  12  0.034550            0.004864
12  13  0.033489            0.004939
13  14  0.033105            0.004799
14  15  0.032704            0.004782
15  16  0.031995            0.004899
16  17  0.031337            0.004575
17  18  0.030109            0.005185
18  19  0.029680            0.005192
19  20  0.029223            0.004962
rms_cca_distance:
     p      mean  standard_deviation
0    1  0.525533            0.256161
1    2  0.069504            0.012882
2    3  0.088734            0.009029
3    4  0.252361            0.081862
4    5  0.126856            0.013084
5    6  0.141515            0.012029
6    7  0.221383            0.055702
7    8  0.180711            0.022460
8    9  0.210041            0.024030
9   10  0.230481            0.026864
10  11  0.244004            0.040005
11  12  0.283121            0.027968
12  13  0.296309            0.027656
13  14  0.330686            0.031417
14  15  0.334938            0.030241
15  16  0.359512            0.021286
16  17  0.357224            0.029239
17  18  0.390234            0.019264
18  19  0.402841            0.019841
19  20  0.404040            0.016546
n_principal_paths_candidates:
[3, 6, 9, 18]

treefit.treefit() returns a treefit.fit.Fit object that summarizes the analysis of Treefit. We will explain how to interpret the results in the next section. For now, we may focus on learning how to use Treefit.

As we will see later, it is helpful to visualize the results using treefit.plot(). By executing treefit.plot(fit_tree_like), we can obtain the following two user-friendly visual plots, which makes it easier to interpret the results of the Treefit analysis.

In [11]: treefit.plot(fit_tree_like)
_images/plot_tree_like_data_estimation_result.png

Figure 3. The output for the tree-like data shown in Figure 1

2.2. Less tree-like, noisy toy data

We can analyze the second toy data in the same manner.

In [12]: fit_less_tree_like = \
   ....:     treefit.treefit({'expression': star_less_tree_like},
   ....:                     name='less-tree-like')
   ....: 

In [13]: print(fit_less_tree_like)
treefit.fit.Fit: less-tree-like
max_cca_distance:
     p      mean  standard_deviation
0    1  0.805926            0.226182
1    2  0.273711            0.061454
2    3  0.231613            0.039771
3    4  0.210874            0.035432
4    5  0.196933            0.032405
5    6  0.182219            0.026901
6    7  0.173168            0.028497
7    8  0.164329            0.028447
8    9  0.150368            0.029421
9   10  0.139997            0.025840
10  11  0.136604            0.026209
11  12  0.122230            0.022835
12  13  0.115212            0.023572
13  14  0.111331            0.021646
14  15  0.108842            0.020588
15  16  0.102107            0.018567
16  17  0.098205            0.019104
17  18  0.095000            0.020071
18  19  0.089316            0.018986
19  20  0.085445            0.014010
rms_cca_distance:
     p      mean  standard_deviation
0    1  0.805926            0.226182
1    2  0.460137            0.047803
2    3  0.577808            0.068314
3    4  0.627684            0.045517
4    5  0.610238            0.041886
5    6  0.606020            0.052572
6    7  0.600699            0.052121
7    8  0.591505            0.045530
8    9  0.569942            0.051620
9   10  0.555651            0.044463
10  11  0.553027            0.050516
11  12  0.550825            0.051761
12  13  0.555187            0.046762
13  14  0.558086            0.038971
14  15  0.568755            0.036608
15  16  0.565789            0.033530
16  17  0.572950            0.027621
17  18  0.567119            0.030153
18  19  0.569511            0.028331
19  20  0.577514            0.026444
n_principal_paths_candidates:
[3, 13, 17, 19]
In [14]: treefit.plot(fit_less_tree_like)
_images/plot_less_tree_like_data_estimation_result.png

Figure 4. The output for the noisy data shown in Figure 2

2.3. Comparing two results

We can compare different results by passing all results to treefit.plot() as follows:

In [15]: treefit.plot(fit_tree_like, fit_less_tree_like)
_images/plot_estimation_results.png

Figure 5. Comparison of the plots shown in Figure 3 and Figure 4

3. Interpreting the results

Before interpreting the previous results, we briefly summarize the process of the Treefit analysis that consists of the following three steps.

  1. First, Treefit repeatedly “perturbs” the input data (i.e., adds some small noise to the original row count data or normalized expression data) in order to produce many slightly different datasets that may have been acquired in the biological experiment.

  2. Second, for each dataset, Treefit calculates a distance matrix that represents the dissimilarities between sample cells and then constructs a tree from each distance matrix. The current version of Treefit computes a minimum spanning tree (MST) that has been widely used for trajectory inference.

  3. Finally, Treefit evaluates the goodness-of-fit between the data and tree models. The underlying idea of this method is that the structure of trees inferred from tree-like data tends to have high robustness to noise, compared to non-tree-like data. Therefore, Treefit measures the mutual similarity between estimated trees in order to check the stability of the tree structures. To this end, Treefit constructs a <i>p</i>-dimensional subspace that extracts the main features of each tree structure and then measuring mutual similarities between the subspaces by using a special type of metrics called the Grassmann distance. In principle, when the estimated trees are mutually similar in their structure, the mean and standard deviation (SD) of the Grassmann distance are small.

Although the word “Grassmann distance” may sound so unfamiliar to some readers, the concept appears in different disguises in various practical contexts. For example, the Grassmann distance has a close connection to canonical correlation analysis (CCA). Treefit provides two Grassmann distances treefit.fit.Fit.max_cca_distance and treefit.fit.Fit.rms_cca_distance that can be used for different purposes as we now explain.

3.1. Estimating the tree-likeness of data

The Treefit analysis using the first Grassmann distance treefit.fit.Fit.max_cca_distance (shown in the left panel of Figure 5) tells us the goodness-of-fit between data and tree models. In principle, as mentioned earlier, if the mean and SD of treefit.fit.Fit.max_cca_distance are small, then this means that the estimated trees are mutually similar in their structure. As can be observed, the distance changes according to the dimensionality p of the feature space, but treefit.fit.Fit.max_cca_distance has the property that the value decreases monotonically as p increases for any datasets.

Comparing the Treefit results for the two datasets, we see that the mean Grassmann distance for the first data does not fall below the second one regardless of the value of p and that the SD of the Grassmann distance for the first data is very small compared to the second data. These results imply that the estimated tree structures are very robust to noise in the first case but not in the second case. Thus, Treefit has verified that the first data are highly tree-like while the second data are not.

3.2. Predicting the number of principal paths in the underlying tree

The Treefit analysis using the other Grassmann distance treefit.fit.Fit.rms_cca_distance (shown in the right panel of Figure 5) is useful to infer the number of “principal paths” in the best-fit tree. From a biological perspective, this analysis can be used to discover a novel or unexpected cell type from single-cell gene expression.

Unlike the previous Grassmann distance, the mean value of treefit.fit.Fit.rms_cca_distance can fluctuate depending on the value of p. Interestingly, we can predict the number of principal paths in the best-fit tree by exploring for which p the distance value reaches “the bottom of a valley” (i.e., attains a local minimum). More precisely, when treefit.fit.Fit.rms_cca_distance attains a local minimum at a certain p, the value p +1 indicates the number of principal paths in the best-fit tree. treefit.fit.Fit.n_principal_paths_candidates has these p +1 values. We don’t need to calculate them manually. When Treefit produces a plot having more than one valleys, the smallest p is usually most informative for the prediction. The smallest p +1 value can be obtained by treefit.fit.Fit.n_principal_paths_candidates[0].

Comparing the Treefit results for the two datasets, we first see that both plots attains a local minimum at p =2. This means that for both datasets the best-fit tree has p +1=3 principle paths, which is correct because both were generated from the same star tree with three arms. Another important point to be made is that the SD of the Grassmann distance for the first data is very small at p =2 compared to that for the second data; in other words, Treefit made this prediction more confidently for the first dataset than for the second one. This result is reasonable because the first dataset is much less noisy than the second one. Thus, Treefit has correctly determined the number of principal paths in the underlying tree together with the goodness-of-fit for each dataset.