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.

The PseudotimeKernel recovers differentiation directionality: We infuse directionality into a k-NN graph using any pseudotime; edges that point into the pseudotemporal “past” are donw-weighted. This can be understood as an adapted, soft version of the Palantir algorithm [Setty et al., 2019].

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:

Removing CLPs does not resolve the problem: Even with CLPs removed, projected velocities still point oppposite to what’s known. The top-likelihood genes look different now - the CLP outliers have been removed; however, many of these top-influential genes show signs of state and time-dependent kinetic parameters. For example, in ANK1, Erythroid cells seem to require their own parameter set and in RPS16, the direction is inverted, i.e., an up-regulation is detected as a down-regulation, probably due to transcriptional bursting [Barile et al., 2021, Bergen et al., 2021].

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