Interactive version

Pancreas Basics

This tutorial shows how to apply CellRank in order to infer initial or terminal states of a developmental process and how to compute probabilistic fate mappings. 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. (2018). For more info on scVelo, see the documentation or read the article.

This is the high level mode to interact with CellRank, which is quick and simple but does not offer as many options as the lower level mode, in which you interact directly with our kernels and estimators. If you would like to get to know this more advanced way of interacting with CellRank, see the pancreas advanced tutorial.

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.

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

First, we need to get the data. The following commands will download the adata object and save it under datasets/endocrinogenesis_day15.5.h5ad.

adata = cr.datasets.pancreas()
scv.utils.show_proportions(adata)
adata

Abundance of ['spliced', 'unspliced']: [0.81 0.19]

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.

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=30, n_neighbors=30)

Filtered out 22024 genes that are detected 20 counts (shared).
Normalized count data: X, spliced, unspliced.
Exctracted 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. The first step, estimating the parameters of the dynamical model, may take a while (~10min). To make sure we only have to run this once, we developed a caching extension called scachepy. scachepy does not only work for recover_dynamics, but it can cache the output of almost any scanpy or scvelo function. To install it, simply run

pip install git+https://github.com/theislab/scachepy

If you don’t want to install scachepy now, don’t worry, the below cell will run without it as well and this is the only place in this tutorial where we’re using it.

try:
    import scachepy
    c = scachepy.Cache('../../cached_files/basic_tutorial/')
    c.tl.recover_dynamics(adata, force=False)
except ModuleNotFoundError:
    print("You don't seem to have scachepy installed, but that's fine, you just have to be a bit patient (~10min). ")
    scv.tl.recover_dynamics(adata)

You don't seem to have scachepy installed, but that's fine, you just have to be a bit patient (~10min).
recovering dynamics
    finished (0:12:15) --> added
    'fit_pars', fitted parameters for splicing dynamics (adata.var)

/opt/hostedtoolcache/Python/3.8.7/x64/lib/python3.8/site-packages/scvelo/tools/dynamical_model.py:683: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  np.array([dm.alpha, dm.beta, dm.gamma, dm.pars[:3]]) / dm.m[-1]
/opt/hostedtoolcache/Python/3.8.7/x64/lib/python3.8/site-packages/scvelo/tools/dynamical_model.py:686: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  np.array([dm.t, dm.tau, dm.t_, dm.pars[4]]) * dm.m[-1]

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.

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
    finished (0:00:03) --> added
    'velocity_graph', sparse matrix with cosine correlations (adata.uns)

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

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

Run CellRank

CellRank builds on the velocities computed by scVelo to compute initial and terminal states of the dynamical process. It further uses the velocities to calculate how likely each cell is to develop from each initial state or towards each terminal state.

Identify terminal states

Terminal states can be computed by running the following command:

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

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 connectivities
    Finish (0:00:00)
Computing eigendecomposition of the transition matrix
Adding `.eigendecomposition`
       `adata.uns['eig_fwd']`
    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 `.eigendecomposition`
       `adata.uns['eig_fwd']`
       `.schur`
       `.schur_matrix`
    Finish (0:00:11)
Computing `3` macrostates
Adding `.macrostates_memberships`
       `.macrostates`
       `.schur`
       `.coarse_T`
       `.coarse_stationary_distribution`
    Finish (0:00:00)
Adding `adata.obs['terminal_states_probs']`
       `adata.obs['terminal_states']`
       `adata.obsm['macrostates_fwd']`
       `.terminal_states_probabilities`
       `.terminal_states`

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”). 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 SLEPSc 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'. 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: additionally to the velocity-based transition probabilities, we use a transition matrix computed on the basis of transcriptomic similarities to make the algorithm more robust. Essentially, we are taking a weighted mean of these two sources of inforamtion, where the weight for transcriptomic similarities is defined by weight_connectivities.

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

cr.pl.terminal_states(adata)

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:

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

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 connectivities
    Finish (0:00:00)
Computing eigendecomposition of the transition matrix
Adding `.eigendecomposition`
       `adata.uns['eig_bwd']`
    Finish (0:00:00)
WARNING: For 1 macrostate, stationary distribution is computed
Adding `.macrostates_memberships`
        `.macrostates`
    Finish (0:00:00)
Adding `adata.obs['initial_states_probs']`
       `adata.obs['initial_states']`
       `adata.obsm['macrostates_bwd']`
       `.terminal_states_probabilities`
       `.terminal_states`

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.

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

Computing lineage probabilities towards terminal states
Computing absorption probabilities
Adding `adata.obsm['to_terminal_states']`
       `adata.obs['to_terminal_states_dp']`
       `.absorption_probabilities`
       `.diff_potential`
    Finish (0:00:00)
Adding lineages to `adata.obsm['to_terminal_states']`
    Finish (0:00:00)

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:

cr.pl.lineages(adata, same_plot=True)

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

Directed PAGA

We can further aggragate the individual fate maps into a cluster-level fate map using an adapted version of PAGA 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.

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:00) --> 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.

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)

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.7/x64/lib/python3.8/site-packages/scvelo/plotting/paga.py:1091: MatplotlibDeprecationWarning:
The inverse_transformed function was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use transformed(transform.inverted()) instead.
  bboxes.append(result.inverse_transformed(ax.transData))
/opt/hostedtoolcache/Python/3.8.7/x64/lib/python3.8/site-packages/scvelo/plotting/paga.py:1091: MatplotlibDeprecationWarning:
The inverse_transformed function was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use transformed(transform.inverted()) instead.
  bboxes.append(result.inverse_transformed(ax.transData))
/opt/hostedtoolcache/Python/3.8.7/x64/lib/python3.8/site-packages/scvelo/plotting/paga.py:1091: MatplotlibDeprecationWarning:
The inverse_transformed function was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use transformed(transform.inverted()) instead.
  bboxes.append(result.inverse_transformed(ax.transData))
/opt/hostedtoolcache/Python/3.8.7/x64/lib/python3.8/site-packages/scvelo/plotting/paga.py:1091: MatplotlibDeprecationWarning:
The inverse_transformed function was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use transformed(transform.inverted()) instead.
  bboxes.append(result.inverse_transformed(ax.transData))
/opt/hostedtoolcache/Python/3.8.7/x64/lib/python3.8/site-packages/scvelo/plotting/paga.py:1091: MatplotlibDeprecationWarning:
The inverse_transformed function was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use transformed(transform.inverted()) instead.
  bboxes.append(result.inverse_transformed(ax.transData))
/opt/hostedtoolcache/Python/3.8.7/x64/lib/python3.8/site-packages/scvelo/plotting/paga.py:1091: MatplotlibDeprecationWarning:
The inverse_transformed function was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use transformed(transform.inverted()) instead.
  bboxes.append(result.inverse_transformed(ax.transData))
/opt/hostedtoolcache/Python/3.8.7/x64/lib/python3.8/site-packages/scvelo/plotting/paga.py:1091: MatplotlibDeprecationWarning:
The inverse_transformed function was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use transformed(transform.inverted()) instead.
  bboxes.append(result.inverse_transformed(ax.transData))

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.

Given the fate maps/probabilistic trajectories, we can ask interesting questions like:

• How does expression of a given gene vary along a specified lineage/trajectory (cellrank.pl.gene_trends and cellrank.pl.heatmap?

• How plastic are different cellular subpopulations (cellrank.pl.cluster_fates(adata, …, mode=’violin’)?

To find out more, check out the CellRank API.

Compute lineage drivers

We can compute the driver genes for all or just the subset of lineages. We can also restric 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.

cr.tl.lineage_drivers(adata)

Computing correlations for lineages `['Epsilon' 'Alpha' 'Beta']` restricted to clusters `None` in layer `X` with `use_raw=False`
Adding `.lineage_drivers`
       `adata.var['to Epsilon corr']`
       `adata.var['to Alpha corr']`
       `adata.var['to Beta corr']`
    Finish (0:00:00)

	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.297290e-07	1.708022e-06	-0.140933	-0.063827	-0.409512	4.725654e-106	6.750934e-104	-0.441431	-0.376558
Anpep	0.456914	7.357066e-136	7.357066e-133	0.425527	0.487202	-0.063882	1.298457e-03	4.477439e-03	-0.102589	-0.024982	-0.228949	1.017801e-31	2.867046e-30	-0.265542	-0.191697
Gm11837	0.449616	6.364047e-131	4.242698e-128	0.417977	0.480167	-0.045986	2.068035e-02	4.854542e-02	-0.084796	-0.007037	-0.238269	2.593631e-34	7.980404e-33	-0.274681	-0.201175
Irx2	0.399584	1.901441e-100	9.507205e-98	0.366325	0.431823	0.517187	3.359556e-182	3.359556e-179	0.488060	0.545163	-0.640866	0.000000e+00	0.000000e+00	-0.663266	-0.617318
Ccnd2	0.384303	3.214070e-92	1.285628e-89	0.350590	0.417020	0.152544	1.074181e-14	1.732551e-13	0.114262	0.190375	-0.351178	6.074690e-76	4.672839e-74	-0.384874	-0.316547
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
Dlk1	-0.323106	1.064602e-63	1.252473e-61	-0.357567	-0.287767	-0.168021	1.477981e-17	2.869865e-16	-0.205637	-0.129910	0.325836	7.865248e-65	4.494427e-63	0.290563	0.360224
Gng12	-0.342606	4.647161e-72	7.149478e-70	-0.376541	-0.307752	-0.330585	7.894530e-67	9.868163e-65	-0.364847	-0.295428	0.462704	7.128450e-140	2.376150e-137	0.431522	0.492782
Nkx6-1	-0.349302	4.411242e-75	7.352071e-73	-0.383051	-0.314622	-0.318417	8.753797e-62	8.753797e-60	-0.352999	-0.282965	0.457423	3.288706e-136	9.396304e-134	0.426055	0.487693
Nnat	-0.357368	7.895266e-79	1.579053e-76	-0.390887	-0.322902	-0.324286	3.462537e-64	4.073572e-62	-0.358716	-0.288976	0.466845	8.497522e-143	3.399009e-140	0.435810	0.496771
Pdx1	-0.370468	3.605169e-85	8.011487e-83	-0.403604	-0.336361	-0.332613	1.076943e-67	1.435924e-65	-0.366821	-0.297507	0.481221	2.648652e-153	1.765768e-150	0.450709	0.510609

2000 rows × 15 columns

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

cr.pl.lineage_drivers(adata, lineage="Alpha", n_genes=5)

Plot smooth gene expression trends along trajectories

The functions demonstrated above are the main functions of CellRank: computing initial and terminal states and probabilistic fate maps. We can use the computed probabilites now to e.g. smooth gene expression trends along lineages.

Let’s start with a temporal ordering for the cells. To get this, we can compute scVelo’s latent time, as computed previously, or alternatively, we can just use CellRank’s initial states to compute a Diffusion Pseudotime.

# compue dpt, starting from CellRank defined root cell
root_idx = np.where(adata.obs['initial_states'] == 'Ngn3 low EP')[0][0]
adata.uns['iroot'] = root_idx
sc.tl.dpt(adata)

scv.pl.scatter(adata, color=['clusters', root_idx, 'latent_time', 'dpt_pseudotime'], fontsize=16,
               cmap='viridis', perc=[2, 98], colorbar=True, rescale_color=[0, 1],
               title=['clusters', 'root cell', 'latent time', 'dpt pseudotime'])

WARNING: Trying to run `tl.dpt` without prior call of `tl.diffmap`. Falling back to `tl.diffmap` with default parameters.

We can plot dynamics of genes in pseudotime along individual trajectories, defined via the fate maps we computed above.

model = cr.ul.models.GAM(adata)
cr.pl.gene_trends(adata, model=model, data_key='X',
                  genes=['Pak3','Neurog3', 'Ghrl'], ncols=3,
                  time_key='latent_time', same_plot=True, hide_cells=True,
                  figsize=(15, 4), n_test_points=200)

Computing trends using `1` core(s)
    Finish (0:00:00)
Plotting trends
/home/runner/work/cellrank_notebooks/cellrank_notebooks/cellrank/cellrank/pl/_utils.py:625: UserWarning: This figure was using constrained_layout==True, but that is incompatible with subplots_adjust and or tight_layout: setting constrained_layout==False.
  fig.tight_layout()

We can also visualize the lineage driver computed above in a heatmap:

cr.pl.heatmap(adata, model, genes=adata.var['to Alpha corr'].sort_values(ascending=False).index[:100],
              show_absorption_probabilities=True,
              lineages="Alpha", n_jobs=1, backend='loky')

Computing trends using `1` core(s)
    Finish (0:00:15)
/home/runner/work/cellrank_notebooks/cellrank_notebooks/cellrank/cellrank/pl/_heatmap.py:457: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  g.fig.add_axes(cax)
/home/runner/work/cellrank_notebooks/cellrank_notebooks/cellrank/cellrank/pl/_heatmap.py:467: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  g.fig.add_axes(cax)

did not converge
did not converge