# CellRank Meets Pseudotime#

## Preliminaries#

In this tutorial, you will learn how to:

compute a pseudotime using

`Diffusion pseudotime (DPT)`

[Haghverdi*et al.*, 2016].use CellRank’s

`PseudotimeKernel`

to compute a transition matrix based on*any*pseudotime of your liking.visualize the transition matrix in a low-dimensional embedding.

Along the way, we’ll see an example where RNA velocity does not work well [Bergen *et al.*, 2020, Bergen *et al.*, 2021, La Manno *et al.*, 2018]; this motivates us to use the `PseudotimeKernel`

.

This tutorial notebook can be downloaded using the following link.

Note

If you want to run this on your own data, you will need:

a scRNA-seq dataset for which you have computed a pseudotime using a tool like

`dpt()`

[Haghverdi*et al.*, 2016], Palantir [Setty*et al.*, 2019] or Slingshot [Street*et al.*, 2018].

Note

If you encounter any bugs in the code, our if you have suggestions for new features, please open an issue. If you have a general question or something you would like to discuss with us, please post on the scverse discourse. You can also contact us using info@cellrank.org.

### Import packages & data#

```
import sys
if "google.colab" in sys.modules:
!pip install -q git+https://github.com/theislab/cellrank
```

```
import numpy as np
import cellrank as cr
import scanpy as sc
import scvelo as scv
scv.settings.verbosity = 3
scv.settings.set_figure_params("scvelo")
sc.settings.set_figure_params(frameon=False, dpi=100)
cr.settings.verbosity = 2
```

```
Global seed set to 0
```

```
import warnings
warnings.simplefilter("ignore", category=UserWarning)
```

To demonstrate the appproach in this tutorial, we will use a scRNA-seq dataset of human bone marrow [Setty *et al.*, 2019], which can be conveniently acessed through `bone_marrow`

.

```
adata = cr.datasets.bone_marrow()
adata
```

```
AnnData object with n_obs × n_vars = 5780 × 27876
obs: 'clusters', 'palantir_pseudotime', 'palantir_diff_potential'
var: 'palantir'
uns: 'clusters_colors', 'palantir_branch_probs_cell_types'
obsm: 'MAGIC_imputed_data', 'X_tsne', 'palantir_branch_probs'
layers: 'spliced', 'unspliced'
```

## Check RNA velocity on this data#

Before diving into the actual `PseudotimeKernel`

, let’s motivate this choice a bit. We’ve seen that RNA velocity works well across a range of datasets including the pancres data from the CellRank Meets RNA Velocity tutorial [Bastidas-Ponce *et al.*, 2019, Bergen *et al.*, 2020]; so let’s check how RNA velocity performs on this dataset.

We’ll check the ratio of spliced to unspliced counts, go through some basic preprocessing, run `scvelo`

, compute a transition matrix using the `VelocityKernel`

and visualize it. To learn more about these steps, please see the CellRank Meets RNA Velocity tutorial.

```
scv.pl.proportions(adata)
```

This looks fine, the percentage of unspliced reads is about what we would expect for 10x Chromium data [La Manno *et al.*, 2018]. Next, filter out genes which don’t have enough spliced/unspliced counts, normalize and log transform the data and restrict to the top highly variable genes. Further, compute principal components and moments for velocity estimation.

```
scv.pp.filter_and_normalize(
adata, min_shared_counts=20, n_top_genes=2000, subset_highly_variable=False
)
sc.tl.pca(adata)
sc.pp.neighbors(adata, n_pcs=30, n_neighbors=30, random_state=0)
scv.pp.moments(adata, n_pcs=None, n_neighbors=None)
```

```
Filtered out 20068 genes that are detected 20 counts (shared).
Normalized count data: X, spliced, unspliced.
Extracted 2000 highly variable genes.
Logarithmized X.
computing moments based on connectivities
finished (0:00:02) --> added
'Ms' and 'Mu', moments of un/spliced abundances (adata.layers)
```

Use the `dynamical model`

from scVelo to estimate model parameters and compute velocities. On my MacBook using 8 cores, the below cell takes about 2 min to execute.

```
scv.tl.recover_dynamics(adata, n_jobs=8)
scv.tl.velocity(adata, mode="dynamical")
```

```
recovering dynamics (using 8/8 cores)
```

```
Global seed set to 0
Global seed set to 0
Global seed set to 0
Global seed set to 0
Global seed set to 0
Global seed set to 0
Global seed set to 0
Global seed set to 0
```

```
finished (0:01:16) --> added
'fit_pars', fitted parameters for splicing dynamics (adata.var)
computing velocities
finished (0:00:01) --> added
'velocity', velocity vectors for each individual cell (adata.layers)
```

Set up the `VelocityKernel`

from the `anndata.AnnData`

object containing the scVelo-computed velocities and compute a cell-cell transition matrix.

```
vk = cr.kernels.VelocityKernel(adata)
vk.compute_transition_matrix()
```

```
Computing transition matrix using `'deterministic'` model
```

```
Using `softmax_scale=1.5903`
```

```
Finish (0:00:02)
```

```
VelocityKernel[n=5780, model='deterministic', similarity='correlation', softmax_scale=1.59]
```

Visualize via stream lines an a t-SNE embedding [Van der Maaten and Hinton, 2008]:

```
vk.plot_projection(basis="tsne")
```

```
Projecting transition matrix onto `tsne`
Adding `adata.obsm['T_fwd_tsne']`
Finish (0:00:00)
```

Note

Arrows point opposite the known differentiation trajectory in which hematopoietic stems cells (HSCs) differentiate via intermediate states towards Monocyotes (Mono), Dendritic cells (DCs), etc [Setty *et al.*, 2019]. That’s not just a result of the low-dimensional representation, feel free to use CellRank to compute initial and terminal states on this data (see Computing Initial and Terminal States tutorial) and you’ll find them to be inconsistend with biological knowledge as well.

To explore why this may be the case, let’s look into the most influential genes driving the velocity flow here:

```
top_genes = adata.var["fit_likelihood"].sort_values(ascending=False).index
scv.pl.scatter(adata, basis=top_genes[:10], ncols=5, frameon=False)
```

In all of the top likelihood genes, the common lymphoid progenitor cells (CLPs) represent an outlier population. Since the current `scvelo`

model does not account for state-dependent kinetic parmeters, this means the CLPs bias the parameter values for all other cells. We explored what happens if we remove CLPs and re-run the above analysis steps:

There’s an easy way in CellRank to overcome these difficulties - use another kernel! In this tutorial, we’ll use the `PseudotimeKernel`

because hematopoiesis is a well-studied system where traditional pseutodime methods work well.

## Use pseudotime to recover directed differentiation#

### Choosing the right pseudotime#

There are many pseudotime algorithms out there, so how do you choose the right one for your data [Saelens *et al.*, 2019]? We’ll do a very superficial analysis here and just compare two methods: `diffusion pseudotime (DPT)`

and the Palantir pseudotime [Haghverdi *et al.*, 2016, Setty *et al.*, 2019].

The Palantir pseudotime has been precomputed for this dataset, check the original tutorial and the `scanpy interface`

to learn how to do this. To compute DPT on this dataset, we’ll start by computing a `diffusion map`

[Coifman *et al.*, 2005, Haghverdi *et al.*, 2015].

```
sc.tl.diffmap(adata)
```

For DPT, we manually have to suply a root cell (recall, we’re not using any RNA velocity here). One (semi-manual) way of doing this is by using extrema of diffusion components:

```
root_ixs = 2394 # has been found using `adata.obsm['X_diffmap'][:, 3].argmax()`
scv.pl.scatter(
adata,
basis="diffmap",
c=["clusters", root_ixs],
legend_loc="right",
components=["2, 3"],
)
adata.uns["iroot"] = root_ixs
```

Once we found a root cell we’re happy with (a cell from the HSC cluster), we can compute DPT and compare it with the precomputed Palantir pseudotime:

```
sc.tl.dpt(adata)
sc.pl.embedding(
adata,
basis="tsne",
color=["dpt_pseudotime", "palantir_pseudotime"],
color_map="gnuplot2",
)
```

```
/Users/marius/miniforge3/envs/py39_arm_cr/lib/python3.9/site-packages/scanpy/plotting/_tools/scatterplots.py:163: MatplotlibDeprecationWarning: The get_cmap function was deprecated in Matplotlib 3.7 and will be removed two minor releases later. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap(obj)`` instead.
cmap = copy(get_cmap(cmap))
```

It seems like DPT is a bit biased towards the CLPs; it assigns very high values to that cluster which masks variation among the other states. We can further explore this with violin plots to visualize the distribution of pseudotime values per cluster, restricted to those clusters we expect to belong to a certain trajectory, e.g. the Monocyte or Erythroid trajectories:

```
mono_trajectory = ["HSC_1", "HSC_2", "Precursors", "Mono_1", "Mono_2"]
ery_trajectory = ["HSC_1", "Ery_1", "Ery_2"]
# plot the Monocyte trajectory
mask = np.in1d(adata.obs["clusters"], mono_trajectory)
sc.pl.violin(
adata[mask],
keys=["dpt_pseudotime", "palantir_pseudotime"],
groupby="clusters",
rotation=-90,
order=mono_trajectory,
)
# plot the Erythroid trajectory
mask = np.in1d(adata.obs["clusters"], ery_trajectory)
sc.pl.violin(
adata[mask],
keys=["dpt_pseudotime", "palantir_pseudotime"],
groupby="clusters",
rotation=-90,
order=ery_trajectory,
)
```

This is a really coarse analysis and only meant to give us a rough idea of which pseudotime to use; generally speaking, they both look good on this dataset. As expected, pseudotimes on average increase as we go towards more mature states. As the Palantir pseudotime appears to be less biased towards CLPs, let’s use that for the `PseudotimeKernel`

below.

### Compute a transition matrix#

Let’s use the Palantir pseudotime to compute a directed cell-cell transition matrix using the `PseudotimeKernel`

:

```
pk = cr.kernels.PseudotimeKernel(adata, time_key="palantir_pseudotime")
pk.compute_transition_matrix()
print(pk)
```

```
Computing transition matrix based on pseudotime`
```

```
Finish (0:00:01)
PseudotimeKernel[n=5780]
```

We can again visualize this transition matrix via streamlines in the t-SNE embedding.

Note

We do not make use of RNA velocity here, CellRank implements a general way of visualizing k-NN graph-based transition matrices via streamlines in any embedding. Thus, the dynamics in the following plot are purely informed by the pseudotime and the k-NN graph, and not by RNA velocity.

```
pk.plot_projection(basis="tsne", recompute=True)
```

```
Projecting transition matrix onto `tsne`
Adding `adata.obsm['T_fwd_tsne']`
Finish (0:00:00)
```

This looks much better, the projected dynamics now agree with what is known from biology.

Note

This is only a low dimensional representation which we shouldn’t trust too much; CellRank contains powerful tools to asses the dynamics in high dimensionional data directly via `estimators`

.

## Closing matters#

### What’s next?#

In this tutorial, you learned how to use CellRank to compute a transition matrix using any precomputed pseudotime and how it can be visualized in low dimensions. The real power of CellRank comes in when you use estimators to analyze the transition matrix directly, rather than projecting it. For the next steps, we recommend to:

go through the initial and terminal states tutorial to learn how to use the transition matrix to automatically identify initial and terminal states.

take a look at the

`API`

to learn about parameter values you can use to adapt these computations to your data.explore the vast amount of pseudotime methods to find the one that works best for your data [Saelens

*et al.*, 2019].

### Package versions#

```
cr.logging.print_versions()
```

```
cellrank==1.5.1+gedbc651e scanpy==1.9.3 anndata==0.8.0 numpy==1.24.2 numba==0.57.0rc1 scipy==1.10.1 pandas==1.5.3 pygpcca==1.0.4 scikit-learn==1.1.3 statsmodels==0.13.5 python-igraph==0.10.4 scvelo==0.3.0 pygam==0.8.0 matplotlib==3.7.0 seaborn==0.12.2
```