Plot lineages

This example shows how to plot the absorption probabilities of the process.

CellRank computes absorption probabilities to estimate how likely each individual cell is to transition into each of the identified terminal or intermediate states. Absorption probabilities refer to the probability of a random walk starting in cell \(i\) to reach terminal state \(j\) before reaching any other terminal state.

Throughout CellRank, we use the terms lineage probabilities, fate probabilities and absorption probabilities interchangeably to describe the same set of probabilities.

import cellrank as cr

adata = cr.datasets.pancreas_preprocessed("../example.h5ad")
adata

Out:

AnnData object with n_obs × n_vars = 2531 × 2000
    obs: 'day', 'proliferation', 'G2M_score', 'S_score', 'phase', 'clusters_coarse', 'clusters', 'clusters_fine', 'louvain_Alpha', 'louvain_Beta', 'initial_size_unspliced', 'initial_size_spliced', 'initial_size', 'n_counts', 'velocity_self_transition', 'dpt_pseudotime'
    var: 'highly_variable_genes', 'gene_count_corr', 'means', 'dispersions', 'dispersions_norm', 'fit_r2', 'fit_alpha', 'fit_beta', 'fit_gamma', 'fit_t_', 'fit_scaling', 'fit_std_u', 'fit_std_s', 'fit_likelihood', 'fit_u0', 'fit_s0', 'fit_pval_steady', 'fit_steady_u', 'fit_steady_s', 'fit_variance', 'fit_alignment_scaling', 'velocity_genes'
    uns: 'clusters_colors', 'clusters_fine_colors', 'diffmap_evals', 'iroot', 'louvain_Alpha_colors', 'louvain_Beta_colors', 'neighbors', 'pca', 'recover_dynamics', 'velocity_graph', 'velocity_graph_neg', 'velocity_params'
    obsm: 'X_diffmap', 'X_pca', 'X_umap', 'velocity_umap'
    varm: 'PCs', 'loss'
    layers: 'Ms', 'Mu', 'fit_t', 'fit_tau', 'fit_tau_', 'spliced', 'unspliced', 'velocity', 'velocity_u'
    obsp: 'connectivities', 'distances'

First, we need to compute the terminal states and the absorption probabilities towards them.

cr.tl.terminal_states(
    adata,
    cluster_key="clusters",
    n_cells=30,
    n_states=3,
    softmax_scale=4,
    show_progress_bar=False,
)
cr.tl.lineages(adata)

Out:

  0%|          | 0/90 [00:00<?, ?/s]
  8%|7         | 7/90 [00:00<00:01, 49.69/s]
 17%|#6        | 15/90 [00:00<00:01, 54.73/s]
 26%|##5       | 23/90 [00:00<00:01, 56.44/s]
 34%|###4      | 31/90 [00:00<00:01, 57.74/s]
 43%|####3     | 39/90 [00:00<00:00, 58.37/s]
 52%|#####2    | 47/90 [00:00<00:00, 58.63/s]
 61%|######1   | 55/90 [00:00<00:00, 58.66/s]
 70%|#######   | 63/90 [00:01<00:00, 56.49/s]
 79%|#######8  | 71/90 [00:01<00:00, 51.60/s]
 87%|########6 | 78/90 [00:01<00:00, 49.22/s]
 94%|#########4| 85/90 [00:01<00:00, 47.47/s]
100%|##########| 90/90 [00:01<00:00, 52.17/s]

We can now plot the terminal states. By default, we plot the degree of membership, which is available only to the cellrank.tl.estimators.GPCCA estimator.

cr.pl.lineages(adata)
absorption probabilities (fwd)

Out:

/home/docs/checkouts/readthedocs.org/user_builds/cellrank/envs/latest/lib/python3.8/site-packages/scvelo/plotting/utils.py:115: MatplotlibDeprecationWarning: Calling gca() with keyword arguments was deprecated in Matplotlib 3.4. Starting two minor releases later, gca() will take no keyword arguments. The gca() function should only be used to get the current axes, or if no axes exist, create new axes with default keyword arguments. To create a new axes with non-default arguments, use plt.axes() or plt.subplot().
  ax = pl.figure(None, figsize, dpi=dpi).gca(projection=projection)

We can also plot only a subset of these lineages or plot the most likely cells.

cr.pl.lineages(adata, ["Alpha", "Beta"], discrete=True)
terminal states (fwd)

Out:

/home/docs/checkouts/readthedocs.org/user_builds/cellrank/envs/latest/lib/python3.8/site-packages/scvelo/plotting/utils.py:115: MatplotlibDeprecationWarning: Calling gca() with keyword arguments was deprecated in Matplotlib 3.4. Starting two minor releases later, gca() will take no keyword arguments. The gca() function should only be used to get the current axes, or if no axes exist, create new axes with default keyword arguments. To create a new axes with non-default arguments, use plt.axes() or plt.subplot().
  ax = pl.figure(None, figsize, dpi=dpi).gca(projection=projection)

Lastly, we can also plot the absorption probabilities separately, one plot for each lineage. By default this also shows the differentiation potential, defined in [Setty19] as the entropy over the absorption probabilities.

cr.pl.lineages(adata, same_plot=False)
to Epsilon, to Alpha, to Beta

To see how to compute and plot the lineage driver genes, refer to Plot potential lineage drivers.

Total running time of the script: ( 0 minutes 24.548 seconds)

Estimated memory usage: 720 MB

Gallery generated by Sphinx-Gallery