Using Velocity and Cellrank with Streams#

This example shows how MIRA’s lineage inference, plotting, and tracing capabilities can be combined with CellRank, which offers a generalized version of the stochastic markov chain model of differentiation that can utilize velocity information. CellRank also provides more advanced and scalable methods for terminal/initial state selection and terminal lineage probability calculations.

This tutorial begins by following the basic CellRank velocity tutorial. After computing absoption probabilities for each cell, we can then use MIRA to parse those probabilities into bifurcating lineages.

Click here to skip to lineage parsing and visualization with MIRA. Otherwise,

ScVelo/CellRank workflow#

[1]:
import scvelo as scv
import scanpy as sc
import cellrank as cr
import numpy as np
import mira
import matplotlib.pyplot as plt

scv.settings.verbosity = 3
cr.settings.verbosity = 2

import warnings
warnings.simplefilter("ignore")
mira.utils.pretty_sderr()
[1]:
[2]:
adata = scv.datasets.pancreas()

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 20801 genes that are detected 20 counts (shared).
Normalized count data: X, spliced, unspliced.
Extracted 2000 highly variable genes.
Logarithmized X.
OMP: Info #271: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
computing moments based on connectivities
    finished (0:00:00) --> added
    'Ms' and 'Mu', moments of un/spliced abundances (adata.layers)
[3]:
scv.tl.recover_dynamics(adata, n_jobs=8)

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

scv.pl.velocity_embedding_stream(
    adata, basis="umap", legend_fontsize=12, title="", smooth=0.8, min_mass=4
)
recovering dynamics (using 8/12 cores)
    finished (0:02:18) --> added
    'fit_pars', fitted parameters for splicing dynamics (adata.var)
computing velocities
    finished (0:00:04) --> added
    'velocity', velocity vectors for each individual cell (adata.layers)
computing velocity graph (using 1/12 cores)
    finished (0:00:09) --> added
    'velocity_graph', sparse matrix with cosine correlations (adata.uns)
computing velocity embedding
    finished (0:00:00) --> added
    'velocity_umap', embedded velocity vectors (adata.obsm)
../_images/notebooks_tutorial_cellrank_3_5.png
[4]:
cr.tl.terminal_states(adata, cluster_key="clusters", weight_connectivities=0,
                     n_states=5, force_recompute=True)
cr.pl.terminal_states(adata)
Computing transition matrix based on logits using `'deterministic'` mode
/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/ipykernel_launcher.py:2: DeprecationWarning: `cellrank.tl.terminal_states` will be removed in version `2.0`. Please use the `cellrank.kernels` or `cellrank.estimators` interface instead.

/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/cellrank/tl/_init_term_states.py:161: DeprecationWarning: `cellrank.tl.transition_matrix` will be removed in version `2.0`. Please use the `cellrank.kernels` or `cellrank.estimators` interface instead.
  **kwargs,
Estimating `softmax_scale` using `'deterministic'` mode
Setting `softmax_scale=5.3689`
    Finish (0:00:05)
Computing Schur decomposition
Mat Object: 1 MPI processes
  type: seqdense
9.9999999999999922e-01 8.0734852408665247e-03 -1.3939720968571084e-02 1.6711722091868716e-02 -2.6729757323936221e-02
0.0000000000000000e+00 9.8776010871620534e-01 -1.0924894189798409e-02 3.0945357751443709e-02 2.8694207884616341e-02
0.0000000000000000e+00 0.0000000000000000e+00 9.8597300904777974e-01 1.2846676832100178e-02 -2.4631237672415930e-03
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.5327128881414147e-01 -1.1276550785823255e-02
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.2593561028091487e-01
Adding `adata.uns['eigendecomposition_fwd']`
       `.schur_vectors`
       `.schur_matrix`
       `.eigendecomposition`
    Finish (0:00:00)
Computing `5` macrostates
INFO:root:Using pre-computed Schur decomposition
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`
../_images/notebooks_tutorial_cellrank_4_9.png
[5]:
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
/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/ipykernel_launcher.py:1: DeprecationWarning: `cellrank.tl.initial_states` will be removed in version `2.0`. Please use the `cellrank.kernels` or `cellrank.estimators` interface instead.
  """Entry point for launching an IPython kernel.
/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/cellrank/tl/_init_term_states.py:161: DeprecationWarning: `cellrank.tl.transition_matrix` will be removed in version `2.0`. Please use the `cellrank.kernels` or `cellrank.estimators` interface instead.
  **kwargs,
Estimating `softmax_scale` using `'deterministic'` mode
Setting `softmax_scale=5.3689`
    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_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`
../_images/notebooks_tutorial_cellrank_5_7.png
[6]:
cr.tl.transition_matrix(adata, weight_connectivities=0.5)
cr.tl.lineages(adata,)
cr.pl.lineages(adata, same_plot=False)
Computing transition matrix based on logits using `'deterministic'` mode
/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/ipykernel_launcher.py:1: DeprecationWarning: `cellrank.tl.transition_matrix` will be removed in version `2.0`. Please use the `cellrank.kernels` or `cellrank.estimators` interface instead.
  """Entry point for launching an IPython kernel.
Estimating `softmax_scale` using `'deterministic'` mode
Setting `softmax_scale=5.3689`
    Finish (0:00:05)
Using a connectivity kernel with weight `0.5`
Computing transition matrix based on `adata.obsp['connectivities']`
    Finish (0:00:00)
Computing absorption probabilities
/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/ipykernel_launcher.py:2: DeprecationWarning: `cellrank.tl.lineages` will be removed in version `2.0`. Please use the `cellrank.kernels` or `cellrank.estimators` interface instead.

Adding `adata.obsm['to_terminal_states']`
       `.absorption_probabilities`
    Finish (0:00:00)
Fatal Python error: GC object already tracked

Current thread 0x000070000ac42000 (most recent call first):
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/traceback.py", line 121 in format_exception
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/traceback.py", line 167 in format_exc
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/multiprocessing/managers.py", line 210 in handle_request
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/threading.py", line 870 in run
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/threading.py", line 926 in _bootstrap_inner
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/threading.py", line 890 in _bootstrap

Thread 0x000070000a73f000 (most recent call first):
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/socket.py", line 212 in accept
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/multiprocessing/connection.py", line 599 in accept
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/multiprocessing/connection.py", line 453 in accept
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/multiprocessing/managers.py", line 178 in accepter
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/threading.py", line 870 in run
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/threading.py", line 926 in _bootstrap_inner
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/threading.py", line 890 in _bootstrap

Thread 0x0000000117f44600 (most recent call first):
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/threading.py", line 300 in wait
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/threading.py", line 552 in wait
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/multiprocessing/managers.py", line 165 in serve_forever
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/multiprocessing/managers.py", line 597 in _run_server
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/multiprocessing/process.py", line 99 in run
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/multiprocessing/process.py", line 297 in _bootstrap
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/multiprocessing/popen_fork.py", line 74 in _launch
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/multiprocessing/popen_fork.py", line 20 in __init__
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/multiprocessing/context.py", line 277 in _Popen
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/multiprocessing/process.py", line 112 in start
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/multiprocessing/managers.py", line 563 in start
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/multiprocessing/context.py", line 56 in Manager
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/cellrank/ul/_parallelize.py", line 101 in wrapper
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/cellrank/tl/_linear_solver.py", line 468 in _solve_lin_system
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/cellrank/tl/estimators/mixins/_absorption_probabilities.py", line 407 in _compute_absorption_probabilities
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/cellrank/tl/estimators/mixins/_absorption_probabilities.py", line 299 in compute_absorption_probabilities
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/cellrank/tl/_lineages.py", line 71 in lineages
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/cellrank/tl/_utils.py", line 1643 in wrapper
  File "/var/folders/26/fjtg4z911nq8pqkhmsbzx88h0000gp/T/ipykernel_50336/3662459673.py", line 2 in <module>
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/IPython/core/interactiveshell.py", line 3524 in run_code
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/IPython/core/interactiveshell.py", line 3444 in run_ast_nodes
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/IPython/core/interactiveshell.py", line 3229 in run_cell_async
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/IPython/core/async_helpers.py", line 78 in _pseudo_sync_runner
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/IPython/core/interactiveshell.py", line 3003 in _run_cell
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/IPython/core/interactiveshell.py", line 2958 in run_cell
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/ipykernel/zmqshell.py", line 533 in run_cell
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/ipykernel/ipkernel.py", line 353 in do_execute
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/ipykernel/kernelbase.py", line 648 in execute_request
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/ipykernel/kernelbase.py", line 353 in dispatch_shell
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/ipykernel/kernelbase.py", line 446 in process_one
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/ipykernel/kernelbase.py", line 457 in dispatch_queue
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/asyncio/events.py", line 88 in _run
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/asyncio/base_events.py", line 1786 in _run_once
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/asyncio/base_events.py", line 541 in run_forever
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/tornado/platform/asyncio.py", line 199 in start
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/ipykernel/kernelapp.py", line 677 in start
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/traitlets/config/application.py", line 846 in launch_instance
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/ipykernel_launcher.py", line 16 in <module>
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/runpy.py", line 85 in _run_code
  File "/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/runpy.py", line 193 in _run_module_as_main
[dsalynchmacport:50432] *** Process received signal ***
[dsalynchmacport:50432] Signal: Abort trap: 6 (6)
[dsalynchmacport:50432] Signal code:  (0)
[dsalynchmacport:50432] [ 0] 0   libsystem_platform.dylib            0x00007ff80af53e2d _sigtramp + 29
[dsalynchmacport:50432] [ 1] 0   ???                                 0x00007fe8081d6dfc 0x0 + 140634545286652
[dsalynchmacport:50432] [ 2] 0   libsystem_c.dylib                   0x00007ff80ae8ad10 abort + 123
[dsalynchmacport:50432] [ 3] 0   python3.7                           0x000000010f7e6e9a fatal_error + 58
[dsalynchmacport:50432] [ 4] 0   python3.7                           0x000000010f7e8103 Py_FatalError + 19
[dsalynchmacport:50432] [ 5] 0   python3.7                           0x000000010f78950d _PyEval_EvalCodeWithName + 5549
[dsalynchmacport:50432] [ 6] 0   python3.7                           0x000000010f65ebdb _PyFunction_FastCallKeywords + 395
[dsalynchmacport:50432] [ 7] 0   python3.7                           0x000000010f796d17 call_function + 183
[dsalynchmacport:50432] [ 8] 0   python3.7                           0x000000010f78f0b7 _PyEval_EvalFrameDefault + 22295
[dsalynchmacport:50432] [ 9] 0   python3.7                           0x000000010f788174 _PyEval_EvalCodeWithName + 532
[dsalynchmacport:50432] [10] 0   python3.7                           0x000000010f65d0ec _PyFunction_FastCallDict + 444
[dsalynchmacport:50432] [11] 0   python3.7                           0x000000010f78f19e _PyEval_EvalFrameDefault + 22526
[dsalynchmacport:50432] [12] 0   python3.7                           0x000000010f65ecf7 _PyFunction_FastCallKeywords + 679
[dsalynchmacport:50432] [13] 0   python3.7                           0x000000010f796d17 call_function + 183
[dsalynchmacport:50432] [14] 0   python3.7                           0x000000010f78f017 _PyEval_EvalFrameDefault + 22135
[dsalynchmacport:50432] [15] 0   python3.7                           0x000000010f65d01a _PyFunction_FastCallDict + 234
[dsalynchmacport:50432] [16] 0   python3.7                           0x000000010f6616ca method_call + 122
[dsalynchmacport:50432] [17] 0   python3.7                           0x000000010f65f46f PyObject_Call + 127
[dsalynchmacport:50432] [18] 0   python3.7                           0x000000010f78f19e _PyEval_EvalFrameDefault + 22526
[dsalynchmacport:50432] [19] 0   python3.7                           0x000000010f65eb32 _PyFunction_FastCallKeywords + 226
[dsalynchmacport:50432] [20] 0   python3.7                           0x000000010f796d17 call_function + 183
[dsalynchmacport:50432] [21] 0   python3.7                           0x000000010f78e7e2 _PyEval_EvalFrameDefault + 20034
[dsalynchmacport:50432] [22] 0   python3.7                           0x000000010f65eb32 _PyFunction_FastCallKeywords + 226
[dsalynchmacport:50432] [23] 0   python3.7                           0x000000010f796d17 call_function + 183
[dsalynchmacport:50432] [24] 0   python3.7                           0x000000010f78e7e2 _PyEval_EvalFrameDefault + 20034
[dsalynchmacport:50432] [25] 0   python3.7                           0x000000010f65d01a _PyFunction_FastCallDict + 234
[dsalynchmacport:50432] [26] 0   python3.7                           0x000000010f6616ca method_call + 122
[dsalynchmacport:50432] [27] 0   python3.7                           0x000000010f65f46f PyObject_Call + 127
[dsalynchmacport:50432] [28] 0   python3.7                           0x000000010f882f5a t_bootstrap + 122
[dsalynchmacport:50432] [29] 0   python3.7                           0x000000010f807ae4 pythread_wrapper + 36
[dsalynchmacport:50432] *** End of error message ***
../_images/notebooks_tutorial_cellrank_6_11.png
[7]:
scv.tl.velocity_pseudotime(adata)
computing terminal states
    identified 1 region of root cells and 1 region of end points .
    finished (0:00:00) --> added
    'root_cells', root cells of Markov diffusion process (adata.obs)
    'end_points', end points of Markov diffusion process (adata.obs)

Lineage parsing with MIRA#

CellRank is essentially a generalization of Palantir that can be used with pseudotime and velocity data. Since MIRA’s core pseudotime trajectory inference functionalities are derived from Palantir, that means MIRA can easily be extended to analyze velocity data!

The cellrank function above computes absorption probabilities for each cell - the probability of diffusing through the markov chain model of the differentiation to each terminal state. MIRA can parse these probabilities to find the bifurcating structure of the data.

In this case, cellrank did not identify the Delta cells as a terminal lineage based on velocity data, so they did not form a branch on the tree.

A notable difference between the annotated clusters (left) and the tree parsing (right) is the tree parsing identified a group of cells, “Epsilon, Alpha”, in orange, that has still not yet committed to either of these fates, as per the velocity vectors.

[8]:
mira.time.get_tree_structure(adata, threshold = 0.712, start_cell = int(adata.obs.initial_states_probs.argmax()),
            pseudotime_key = 'velocity_pseudotime', cellrank = True)

sc.pl.umap(adata, color = ['clusters','tree_states'], palette = 'Set2_r', frameon = False,
          add_outline=False, outline_width=(0.1, 0), outline_color=('lightgrey','white'), size = 50)
INFO:mira.adata_interface.pseudotime:Added key to obs: tree_states
INFO:mira.adata_interface.pseudotime:Added key to uns: tree_state_names
INFO:mira.adata_interface.pseudotime:Added key to uns: connectivities_tree
... storing 'tree_states' as categorical
../_images/notebooks_tutorial_cellrank_9_1.png

Finding the tree structure of the data enables us to explore regulatory dynamics using MIRA’s streamgraphs. Below, we show the expression of marker genes representing transitions in cell identity.

[9]:
cell_mask = adata.obs.tree_states.str.contains('Epsilon|Alpha|Beta') & (adata.obs.velocity_pseudotime > 0.1)
mira.pl.plot_stream(adata[cell_mask],
                    pseudotime_key='velocity_pseudotime',
                   data = ['Pak3','Neurog3','Ghrl'], log_pseudotime=False, palette = 'Set2',
                   scale_features=True, figsize=(7,4), max_bar_height=0.99,
                   clip = 1, min_pseudotime = 0.2, linewidth = 0.5)
plt.show()
../_images/notebooks_tutorial_cellrank_11_0.png

MIRA swarmgraph is useful for plotting discrete values, such as cluster labels. The tree structure of the data largely recapitulates the cluster annotations.

[10]:
hue_order = adata.obs.clusters.cat.categories.values
[11]:
adata.obs['clusters'] = adata.obs.clusters.astype(str)
mira.pl.plot_stream(adata[cell_mask], pseudotime_key='velocity_pseudotime', style = 'swarm', max_swarm_density=500,
                   data = 'clusters', log_pseudotime=False, size = 5, max_bar_height=0.9,
                   hue_order = hue_order,
                   palette = adata.uns['clusters_colors'],
                   scale_features=True, figsize=(7,4), min_pseudotime = 0.2, linecolor = 'lightgrey')
plt.show()
../_images/notebooks_tutorial_cellrank_14_0.png

Streams may also demonstrate RNA velocity. This streagraph shows relative frequencies of unspliced vs. spliced transcripts for gene Cpe:

[12]:
mira.pl.plot_stream(adata[cell_mask], pseudotime_key='velocity_pseudotime', layers=['Ms','Mu'], style = 'line', size = 5,
                   data = ['Cpe','Cpe'], log_pseudotime=False, palette = ['red','black'], window_size=301,
                   scale_features=True, figsize=(7,4), max_bar_height=0.7, clip = 1, min_pseudotime = 0.2)
plt.show()
../_images/notebooks_tutorial_cellrank_16_0.png

Finally, RNA velocity may be combined with MIRA’s trace_differentiation functionality to study routes through which cells diffuse through the differentiation manifold.

[13]:
cr.tl.transition_matrix(adata, weight_connectivities=0) # only use velocity transitions
Computing transition matrix based on logits using `'deterministic'` mode
/Users/alynch/opt/miniconda3/envs/mirarep/lib/python3.7/site-packages/ipykernel_launcher.py:1: DeprecationWarning: `cellrank.tl.transition_matrix` will be removed in version `2.0`. Please use the `cellrank.kernels` or `cellrank.estimators` interface instead.
  """Entry point for launching an IPython kernel.
Estimating `softmax_scale` using `'deterministic'` mode
Setting `softmax_scale=5.3689`
    Finish (0:00:06)
[13]:
<VelocityKernel>

Traces can be used to find ancestral cell populations starting from lineage termini. RNA velocity indicates Epsilon cells follow a less direct path through the differentiation manifold than one might expect:

[14]:
mira.time.trace_differentiation(adata, start_cells= (adata.obs.terminal_states == 'Epsilon').values,
                               direction='forward', transport_map_key='T_bwd', save_name='data/velocity_test.gif',
                               log_prob = True, num_steps=175, steps_per_frame=1, sqrt_time=True, palette='cividis',
                               vmax_quantile=0.95, figsize = (5,3), size = 1, fps=24, add_outline = False,
                               num_preview_frames=3)

mira.utils.show_gif('data/velocity_test.gif')
INFO:mira.pseudotime.backtrace:Creating transport map ...
INFO:mira.pseudotime.backtrace:Tracing ancestral populations ...
../_images/notebooks_tutorial_cellrank_20_1.png
INFO:mira.pseudotime.backtrace:Creating animation ...
INFO:mira.pseudotime.backtrace:Saving animation ...
INFO:matplotlib.animation:Animation.save using <class 'matplotlib.animation.PillowWriter'>
../_images/notebooks_tutorial_cellrank_20_3.png

And some Beta cells appear to descend from Alpha-like cells, branching later in the differentiation:

[15]:
mira.time.trace_differentiation(adata, start_cells= (adata.obs.terminal_states == 'Beta').values,
                               direction='forward', transport_map_key='T_bwd', save_name='data/velocity_test.gif',
                               log_prob = True, num_steps=150, steps_per_frame=1, sqrt_time=True, palette='cividis',
                               vmax_quantile=0.95, figsize = (5,3), size = 1, fps=24, add_outline = False,
                               num_preview_frames=3)

mira.utils.show_gif('data/velocity_test.gif')
INFO:mira.pseudotime.backtrace:Creating transport map ...
INFO:mira.pseudotime.backtrace:Tracing ancestral populations ...
../_images/notebooks_tutorial_cellrank_22_1.png
INFO:mira.pseudotime.backtrace:Creating animation ...
INFO:mira.pseudotime.backtrace:Saving animation ...
INFO:matplotlib.animation:Animation.save using <class 'matplotlib.animation.PillowWriter'>
../_images/notebooks_tutorial_cellrank_22_3.png

Tracing the differentiation forward, one can see cycling states in Ductal cells:

[16]:
mira.time.trace_differentiation(adata, start_cells= [int(adata.obs.initial_states_probs.argmax())],
                               transport_map_key='T_fwd', save_name='data/velocity_test.gif',
                               log_prob = True, num_steps=30, steps_per_frame=1, sqrt_time=False, palette='cividis',
                               vmax_quantile=0.95, figsize = (5,3), size = 1, fps=8, add_outline = False,
                               num_preview_frames=3)

mira.utils.show_gif('data/velocity_test.gif')
INFO:mira.pseudotime.backtrace:Creating transport map ...
INFO:mira.pseudotime.backtrace:Tracing ancestral populations ...
../_images/notebooks_tutorial_cellrank_24_1.png
INFO:mira.pseudotime.backtrace:Creating animation ...
INFO:mira.pseudotime.backtrace:Saving animation ...
INFO:matplotlib.animation:Animation.save using <class 'matplotlib.animation.PillowWriter'>
../_images/notebooks_tutorial_cellrank_24_3.png

CellRank intends to update their infrastructure with new “kernels” - essentially methods with which to assess cell-cell similarity and transition probabilities. Extension of their methods to include exciting new ideas in optimal transport mean MIRA can leverage these new methods de facto.