CellRank basics#

This tutorial introduces you to CellRank’s high level API for computing initial & terminal states and fate probabilities. Once we have the fate probabilities, this tutorial shows you how to use them to plot a directed [Wolf et al., 2019], to compute putative lineage drivers and to visualize smooth gene expression trends. If you want a bit more control over how initial & terminal states and fate probabilities are computed, then you should check out CellRank’s low level API, composed of kernels and estimators. This really isn’t any more complicated than using scikit-learn, so please do check out the Kernels and estimators tutorial.

In this tutorial, we will use RNA velocity and transcriptomic similarity to estimate cell-cell transition probabilities. Using kernels and estimators, you can apply CellRank even without RNA velocity information, check out our CellRank beyond RNA velocity tutorial. CellRank generalizes beyond RNA velocity and is a widely applicable framework to model single-cell data based on the powerful concept of Markov chains.

The first part of this tutorial is very similar to scVelo’s tutorial on pancreatic endocrinogenesis. The data we use here comes from [Bastidas-Ponce et al., 2019]. For more info on scVelo, see the documentation or take a look at [Bergen et al., 2020].

This tutorial notebook can be downloaded using the following link.

Import packages & data#

Easiest way to start is to download Miniconda3 along with the environment file found here. To create the environment, run conda create -f environment.yml.

[1]:

import sys

if "google.colab" in sys.modules:
    !pip install -q git+https://github.com/theislab/cellrank@dev
    !pip install python-igraph

[2]:

import scvelo as scv
import scanpy as sc
import cellrank as cr
import numpy as np

scv.settings.verbosity = 3
scv.settings.set_figure_params("scvelo")
cr.settings.verbosity = 2

[3]:

import warnings

warnings.simplefilter("ignore", category=UserWarning)
warnings.simplefilter("ignore", category=FutureWarning)
warnings.simplefilter("ignore", category=DeprecationWarning)

First, we need to get the data. The following commands will download the adata object and save it under datasets/endocrinogenesis_day15.5.h5ad. We’ll also show the fraction of spliced/unspliced reads, which we need to estimate RNA velocity.

[4]:

adata = cr.datasets.pancreas()
scv.pl.proportions(adata)
adata

100%|██████████| 33.5M/33.5M [00:01<00:00, 27.4MB/s]
../../_images/notebooks_tutorials_cellrank_basics_8_1.png 
[4]:

AnnData object with n_obs × n_vars = 2531 × 27998
    obs: 'day', 'proliferation', 'G2M_score', 'S_score', 'phase', 'clusters_coarse', 'clusters', 'clusters_fine', 'louvain_Alpha', 'louvain_Beta', 'palantir_pseudotime'
    var: 'highly_variable_genes'
    uns: 'clusters_colors', 'clusters_fine_colors', 'day_colors', 'louvain_Alpha_colors', 'louvain_Beta_colors', 'neighbors', 'pca'
    obsm: 'X_pca', 'X_umap'
    layers: 'spliced', 'unspliced'
    obsp: 'connectivities', 'distances'

Pre-process the data#

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. These are standard scanpy/scvelo functions, for more information about them, see the scVelo API.

[5]:

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

Filtered out 22024 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:00) --> added
    'Ms' and 'Mu', moments of un/spliced abundances (adata.layers)

Run scVelo#

We will use the dynamical model from scVelo to estimate the velocities. Please make sure to have at least version 0.2.3 of scVelo installed to make use parallelisation in scv.tl.recover_dynamics. On my laptop, using 8 cores, the below cell takes about 1:30 min to execute.

[6]:

scv.tl.recover_dynamics(adata, n_jobs=8)

recovering dynamics (using 2/2 cores)
WARNING: Unable to create progress bar. Consider installing `tqdm` as `pip install tqdm` and `ipywidgets` as `pip install ipywidgets`,
or disable the progress bar using `show_progress_bar=False`.
    finished (0:03:07) --> added
    'fit_pars', fitted parameters for splicing dynamics (adata.var)

Once we have the parameters, we can use these to compute the velocities and the velocity graph. The velocity graph is a weighted graph that specifies how likely two cells are to transition into another, given their velocity vectors and relative positions.

[7]:

scv.tl.velocity(adata, mode="dynamical")
scv.tl.velocity_graph(adata)

computing velocities
    finished (0:00:02) --> added
    'velocity', velocity vectors for each individual cell (adata.layers)
computing velocity graph (using 1/2 cores)
    finished (0:00:04) --> added
    'velocity_graph', sparse matrix with cosine correlations (adata.uns)

[8]:

scv.pl.velocity_embedding_stream(
    adata, basis="umap", legend_fontsize=12, title="", smooth=0.8, min_mass=4
)

computing velocity embedding
    finished (0:00:00) --> added
    'velocity_umap', embedded velocity vectors (adata.obsm)
../../_images/notebooks_tutorials_cellrank_basics_17_1.png

Run CellRank#

CellRank offers various ways to infuse directionality into single-cell data. Here, the directional information comes from RNA velocity, and we use this information to compute initial & terminal states as well as fate probabilities for the dynamical process of pancreatic development.

Identify terminal states#

Terminal states can be computed by running the following command:

[9]:

cr.tl.terminal_states(adata, cluster_key="clusters", weight_connectivities=0.2)

Accessing `adata.obsp['T_fwd']`
Computing transition matrix based on logits using `'deterministic'` mode
Estimating `softmax_scale` using `'deterministic'` mode
Setting `softmax_scale=3.7951`
    Finish (0:00:10)
Using a connectivity kernel with weight `0.2`
Computing transition matrix based on `adata.obsp['connectivities']`
    Finish (0:00:00)
Computing eigendecomposition of the transition matrix
Adding `adata.uns['eigendecomposition_fwd']`
       `.eigendecomposition`
    Finish (0:00:00)
WARNING: Unable to import `petsc4py` or `slepc4py`. Using `method='brandts'`
WARNING: For `method='brandts'`, dense matrix is required. Densifying
Computing Schur decomposition
Adding `adata.uns['eigendecomposition_fwd']`
       `.schur_vectors`
       `.schur_matrix`
       `.eigendecomposition`
    Finish (0:00:06)
Computing `3` macrostates
Adding `.macrostates`
       `.macrostates_memberships`
       `.coarse_T`
       `.coarse_initial_distribution
       `.coarse_stationary_distribution`
       `.schur_vectors`
       `.schur_matrix`
       `.eigendecomposition`
    Finish (0:00:00)
Adding `adata.obs['terminal_states']`
       `adata.obs['terminal_states_probs']`
       `.terminal_states`
       `.terminal_states_probabilities`
       `.terminal_states_memberships
    Finish`

/tmp/ipykernel_1982/1878820466.py:1: DeprecationWarning: `cellrank.tl.terminal_states` will be removed in version `2.0`. Please use the `cellrank.kernels` or `cellrank.estimators` interface instead.
  cr.tl.terminal_states(adata, cluster_key="clusters", weight_connectivities=0.2)
/home/runner/work/cellrank_notebooks/cellrank_notebooks/cellrank/cellrank/tl/_init_term_states.py:156: DeprecationWarning: `cellrank.tl.transition_matrix` will be removed in version `2.0`. Please use the `cellrank.kernels` or `cellrank.estimators` interface instead.
  kernel = transition_matrix(
100%|██████████| 2531/2531 [00:08<00:00, 282.11cell/s]
100%|██████████| 2531/2531 [00:01<00:00, 1539.40cell/s]

The most important parameters in the above function are:

  • estimator: this determines what’s going to behind the scenes to compute the terminal states. Options are cr.tl.estimators.CFLARE (“Clustering and Filtering of Left and Right Eigenvectors”) or cr.tl.estimators.GPCCA (“Generalized Perron Cluster Cluster Analysis, [Reuter et al., 2018] and [Reuter et al., 2019], see also our pyGPCCA implementation). The latter is the default, it computes terminal states by coarse graining the velocity-derived Markov chain into a set of macrostates that represent the slow-time scale dynamics of the process, i.e. it finds the states that you are unlikely to leave again, once you have entered them.

  • cluster_key: takes a key from adata.obs to retrieve pre-computed cluster labels, i.e. ‘clusters’ or ‘louvain’. These labels are then mapped onto the set of terminal states, to associate a name and a color with each state.

  • n_states: number of expected terminal states. This parameter is optional - if it’s not provided, this number is estimated from the so-called ‘eigengap heuristic’ of the spectrum of the transition matrix.

  • method: This is only relevant for the estimator GPCCA. It determines the way in which we compute and sort the real Schur decomposition. The default, krylov, is an iterative procedure that works with sparse matrices which allows the method to scale to very large cell numbers. It relies on the libraries SLEPc and PETSc, which you will have to install separately, see our installation instructions. If your dataset is small (<5k cells), and you don’t want to install these at the moment, use method='brandts' [Brandts, 2002]. The results will be the same, the difference is that brandts works with dense matrices and won’t scale to very large cells numbers.

  • weight_connectivities: weight given to cell-cell similarities to account for noise in velocity vectors.

When running the above command, CellRank adds a key terminal_states to adata.obs and the result can be plotted as:

[10]:

cr.pl.terminal_states(adata)
../../_images/notebooks_tutorials_cellrank_basics_24_0.png

Identify initial states#

The same sort of analysis can now be repeated for the initial states, only that we use the function cr.tl.initial_states this time:

[11]:

cr.tl.initial_states(adata, cluster_key="clusters")
cr.pl.initial_states(adata, discrete=True)

Accessing `adata.obsp['T_bwd']`
Computing transition matrix based on logits using `'deterministic'` mode
Estimating `softmax_scale` using `'deterministic'` mode
Setting `softmax_scale=3.7951`
    Finish (0:00:07)
Using a connectivity kernel with weight `0.2`
Computing transition matrix based on `adata.obsp['connectivities']`
    Finish (0:00:00)
Computing eigendecomposition of the transition matrix
Adding `adata.uns['eigendecomposition_bwd']`
       `.eigendecomposition`
    Finish (0:00:00)
For 1 macrostate, stationary distribution is computed
Adding `.macrostates`
       `.macrostates_memberships`
       `.coarse_T`
       `.coarse_initial_distribution
       `.coarse_stationary_distribution`
       `.schur_vectors`
       `.schur_matrix`
       `.eigendecomposition`
    Finish (0:00:00)
Adding `adata.obs['initial_states']`
       `adata.obs['initial_states_probs']`
       `.terminal_states`
       `.terminal_states_probabilities`
       `.terminal_states_memberships
    Finish`

/tmp/ipykernel_1982/1187097337.py:1: DeprecationWarning: `cellrank.tl.initial_states` will be removed in version `2.0`. Please use the `cellrank.kernels` or `cellrank.estimators` interface instead.
  cr.tl.initial_states(adata, cluster_key="clusters")
/home/runner/work/cellrank_notebooks/cellrank_notebooks/cellrank/cellrank/tl/_init_term_states.py:156: DeprecationWarning: `cellrank.tl.transition_matrix` will be removed in version `2.0`. Please use the `cellrank.kernels` or `cellrank.estimators` interface instead.
  kernel = transition_matrix(
100%|██████████| 2531/2531 [00:04<00:00, 545.09cell/s]
100%|██████████| 2531/2531 [00:02<00:00, 918.25cell/s]
../../_images/notebooks_tutorials_cellrank_basics_27_2.png

We found one initial state, located in the Ngn3 low EP cluster.

Compute fate maps#

Once we know the terminal states, we can compute associated fate maps - for each cell, we ask how likely is the cell to develop towards each of the identified terminal states.

[12]:

cr.tl.lineages(adata)
cr.pl.lineages(adata, same_plot=False)

Computing absorption probabilities
WARNING: Unable to import petsc4py. For installation, please refer to: https://petsc4py.readthedocs.io/en/stable/install.html.
Defaulting to `'gmres'` solver.
Adding `adata.obsm['to_terminal_states']`
       `.absorption_probabilities`
    Finish (0:00:00)

/tmp/ipykernel_1982/1878844763.py:1: DeprecationWarning: `cellrank.tl.lineages` will be removed in version `2.0`. Please use the `cellrank.kernels` or `cellrank.estimators` interface instead.
  cr.tl.lineages(adata)
100%|██████████| 3/3 [00:00<00:00, 34.12/s]
../../_images/notebooks_tutorials_cellrank_basics_31_2.png

We can aggregate the above into a single, global fate map where we associate each terminal state with color and use the intensity of that color to show the fate of each individual cell:

[13]:

cr.pl.lineages(adata, same_plot=True)
../../_images/notebooks_tutorials_cellrank_basics_33_0.png

This shows that the dominant terminal state at E15.5 is Beta, consistent with known biology, see e.g. [Bastidas-Ponce et al., 2019].

Directed PAGA#

We can further aggregate the individual fate maps into a cluster-level fate map using an adapted version of [Wolf et al., 2019] with directed edges. We first compute scVelo’s latent time with CellRank identified root_key and end_key, which are the probabilities of being an initial or a terminal state, respectively.

[14]:

scv.tl.recover_latent_time(
    adata, root_key="initial_states_probs", end_key="terminal_states_probs"
)

computing latent time using initial_states_probs, terminal_states_probs as prior
    finished (0:00:01) --> added
    'latent_time', shared time (adata.obs)

Next, we can use the inferred pseudotime along with the initial and terminal states probabilities to compute the directed PAGA.

[15]:

scv.tl.paga(
    adata,
    groups="clusters",
    root_key="initial_states_probs",
    end_key="terminal_states_probs",
    use_time_prior="velocity_pseudotime",
)

running PAGA using priors: ['velocity_pseudotime', 'initial_states_probs', 'terminal_states_probs']
    finished (0:00:00) --> added
    'paga/connectivities', connectivities adjacency (adata.uns)
    'paga/connectivities_tree', connectivities subtree (adata.uns)
    'paga/transitions_confidence', velocity transitions (adata.uns)

[16]:

cr.pl.cluster_fates(
    adata,
    mode="paga_pie",
    cluster_key="clusters",
    basis="umap",
    legend_kwargs={"loc": "top right out"},
    legend_loc="top left out",
    node_size_scale=5,
    edge_width_scale=1,
    max_edge_width=4,
    title="directed PAGA",
)

/opt/hostedtoolcache/Python/3.8.12/x64/lib/python3.8/site-packages/networkx/drawing/nx_pylab.py:717: MatplotlibDeprecationWarning: Passing *transOffset* without *offsets* has no effect. This behavior is deprecated since 3.5 and in 3.6, *transOffset* will begin having an effect regardless of *offsets*. In the meantime, if you wish to set *transOffset*, call collection.set_offset_transform(transOffset) explicitly.
  edge_collection = mpl.collections.LineCollection(
../../_images/notebooks_tutorials_cellrank_basics_40_1.png

We use pie charts to show cell fates averaged per cluster. Edges between clusters are given by transcriptomic similarity between the clusters, just as in normal PAGA.

Compute lineage drivers#

We can compute the driver genes for all or just a subset of lineages. We can also restrict this to some subset of clusters by specifying clusters=... (not shown below). In the resulting dataframe, we also see the p-value, the corrected p-value (q-value) and the 95% confidence interval for the correlation statistic.

[17]:

cr.tl.lineage_drivers(adata)

Adding `adata.varm['terminal_lineage_drivers']`
       `.lineage_drivers`
    Finish (0:00:00)

/tmp/ipykernel_1982/1721556810.py:1: DeprecationWarning: `cellrank.tl.lineage_drivers` will be removed in version `2.0`. Please use the `cellrank.kernels` or `cellrank.estimators` interface instead.
  cr.tl.lineage_drivers(adata)

[17]:
Epsilon_corr Epsilon_pval Epsilon_qval Epsilon_ci_low Epsilon_ci_high Alpha_corr Alpha_pval Alpha_qval Alpha_ci_low Alpha_ci_high Beta_corr Beta_pval Beta_qval Beta_ci_low Beta_ci_high
index
Ghrl 0.802445 0.000000e+00 0.000000e+00 0.788123 0.815898 -0.102534 2.297271e-07 1.708008e-06 -0.140933 -0.063827 -0.409512 4.724494e-106 6.749276e-104 -0.441431 -0.376559
Anpep 0.456913 7.362724e-136 7.362724e-133 0.425527 0.487202 -0.063882 1.298457e-03 4.477439e-03 -0.102589 -0.024982 -0.228948 1.018126e-31 2.867961e-30 -0.265541 -0.191696
Gm11837 0.449617 6.363291e-131 4.242194e-128 0.417977 0.480167 -0.045986 2.068030e-02 4.854529e-02 -0.084796 -0.007037 -0.238269 2.593336e-34 7.979495e-33 -0.274681 -0.201175
Irx2 0.399584 1.900072e-100 9.500358e-98 0.366326 0.431823 0.517187 3.359790e-182 3.359790e-179 0.488060 0.545163 -0.640866 0.000000e+00 0.000000e+00 -0.663266 -0.617318
Ccnd2 0.384302 3.215060e-92 1.286024e-89 0.350590 0.417020 0.152544 1.074184e-14 1.732555e-13 0.114262 0.190375 -0.351177 6.078981e-76 4.676139e-74 -0.384873 -0.316546
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
Dlk1 -0.323105 1.065363e-63 1.253368e-61 -0.357566 -0.287766 -0.168021 1.477933e-17 2.869774e-16 -0.205637 -0.129910 0.325835 7.866581e-65 4.495189e-63 0.290562 0.360224
Gng12 -0.342605 4.652637e-72 7.157903e-70 -0.376540 -0.307751 -0.330585 7.893886e-67 9.867357e-65 -0.364847 -0.295428 0.462703 7.138827e-140 2.379609e-137 0.431521 0.492781
Nkx6-1 -0.349302 4.412816e-75 7.354693e-73 -0.383050 -0.314622 -0.318417 8.753770e-62 8.753770e-60 -0.352999 -0.282965 0.457422 3.293317e-136 9.409478e-134 0.426054 0.487693
Nnat -0.357367 7.903873e-79 1.580775e-76 -0.390886 -0.322901 -0.324287 3.462263e-64 4.073251e-62 -0.358716 -0.288976 0.466844 8.508331e-143 3.403332e-140 0.435809 0.496770
Pdx1 -0.370467 3.608772e-85 8.019494e-83 -0.403603 -0.336360 -0.332613 1.076877e-67 1.435836e-65 -0.366821 -0.297507 0.481220 2.652234e-153 1.768156e-150 0.450709 0.510608

2000 rows × 15 columns

Afterwards, we can plot the top 5 driver genes (based on the correlation), e.g. for the Alpha lineage:

[18]:

cr.pl.lineage_drivers(adata, lineage="Alpha", n_genes=5)
../../_images/notebooks_tutorials_cellrank_basics_46_0.png

What’s next?#

Congratulations! You have successfully gone through some first computations with CellRank. If you want to learn more, you can check out:

  • our low level API, unlocking the full potential of CellRank trough kernels and estimators, flexible classes that compute cell-cell transition probabilities (kernels) and aggregate these to formulate hypothesis about the underlying dynamics (estimators). See the Kernels and estimators tutorial.

  • how CellRank can be used without RNA velocity. See the CellRank beyond RNA velocity tutorial.

  • CellRank external, our interface to third-party libraries, giving you even more possibilities to model single-cell data based on Markov chains, conveniently though the CellRank interface.