Plot projection

This example shows how to project a transition matrix onto a low-dimensional embedding.

import scvelo as scv
import cellrank as cr

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

Out:

  0%|          | 0.00/140M [00:00<?, ?B/s]
  0%|          | 56.0k/140M [00:00<06:22, 384kB/s]
  0%|          | 176k/140M [00:00<03:50, 636kB/s]
  0%|          | 680k/140M [00:00<01:16, 1.90MB/s]
  2%|1         | 2.62M/140M [00:00<00:22, 6.50MB/s]
  4%|4         | 5.65M/140M [00:00<00:12, 11.7MB/s]
  6%|6         | 8.68M/140M [00:00<00:09, 14.9MB/s]
  8%|8         | 11.7M/140M [00:01<00:08, 16.8MB/s]
 11%|#         | 14.7M/140M [00:01<00:07, 18.1MB/s]
 13%|#2        | 17.8M/140M [00:01<00:06, 19.1MB/s]
 15%|#4        | 20.8M/140M [00:01<00:06, 19.5MB/s]
 17%|#7        | 23.9M/140M [00:01<00:06, 19.9MB/s]
 19%|#9        | 26.9M/140M [00:01<00:05, 20.2MB/s]
 21%|##1       | 29.9M/140M [00:01<00:05, 20.4MB/s]
 24%|##3       | 32.9M/140M [00:02<00:05, 20.6MB/s]
 26%|##5       | 36.0M/140M [00:02<00:05, 20.7MB/s]
 28%|##7       | 39.0M/140M [00:02<00:05, 20.7MB/s]
 30%|##9       | 42.0M/140M [00:02<00:04, 20.8MB/s]
 32%|###2      | 45.1M/140M [00:02<00:04, 21.1MB/s]
 34%|###4      | 48.1M/140M [00:02<00:04, 21.1MB/s]
 36%|###6      | 51.1M/140M [00:03<00:04, 20.7MB/s]
 39%|###8      | 54.2M/140M [00:03<00:04, 20.7MB/s]
 41%|####      | 57.2M/140M [00:03<00:04, 20.8MB/s]
 43%|####2     | 60.2M/140M [00:03<00:04, 20.8MB/s]
 45%|####5     | 63.3M/140M [00:03<00:03, 20.8MB/s]
 47%|####7     | 66.3M/140M [00:03<00:03, 20.8MB/s]
 49%|####9     | 69.3M/140M [00:03<00:03, 20.9MB/s]
 52%|#####1    | 72.4M/140M [00:04<00:03, 21.4MB/s]
 54%|#####3    | 75.4M/140M [00:04<00:03, 21.3MB/s]
 56%|#####5    | 78.4M/140M [00:04<00:03, 20.6MB/s]
 58%|#####8    | 81.5M/140M [00:04<00:02, 20.7MB/s]
 60%|######    | 84.5M/140M [00:04<00:02, 20.7MB/s]
 62%|######2   | 87.5M/140M [00:04<00:02, 20.8MB/s]
 65%|######4   | 90.6M/140M [00:05<00:02, 20.8MB/s]
 67%|######6   | 93.6M/140M [00:05<00:02, 20.9MB/s]
 69%|######8   | 96.6M/140M [00:05<00:02, 20.9MB/s]
 71%|#######1  | 99.7M/140M [00:05<00:02, 21.1MB/s]
 73%|#######3  | 103M/140M [00:05<00:01, 21.1MB/s]
 75%|#######5  | 106M/140M [00:05<00:01, 20.7MB/s]
 78%|#######7  | 109M/140M [00:05<00:01, 20.8MB/s]
 80%|#######9  | 112M/140M [00:06<00:01, 20.7MB/s]
 82%|########1 | 115M/140M [00:06<00:01, 20.8MB/s]
 84%|########4 | 118M/140M [00:06<00:01, 20.8MB/s]
 86%|########6 | 121M/140M [00:06<00:00, 20.8MB/s]
 88%|########8 | 124M/140M [00:06<00:00, 20.9MB/s]
 91%|######### | 127M/140M [00:06<00:00, 21.0MB/s]
 93%|#########2| 130M/140M [00:06<00:00, 21.0MB/s]
 95%|#########4| 133M/140M [00:07<00:00, 20.8MB/s]
 97%|#########7| 136M/140M [00:07<00:00, 20.8MB/s]
 99%|#########9| 139M/140M [00:07<00:00, 20.9MB/s]
100%|##########| 140M/140M [00:07<00:00, 19.7MB/s]

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 create the kernel using the high-level function cellrank.tl.transition_matrix(), but any instance of cellrank.tl.kernels.Kernel would do.

k = cr.tl.transition_matrix(
    adata, weight_connectivities=0.2, softmax_scale=4, show_progress_bar=False
)
k

Out:

/home/docs/checkouts/readthedocs.org/user_builds/cellrank/checkouts/latest/examples/kernels/plot_projection.py:17: DeprecationWarning: `cellrank.tl.transition_matrix` will be removed in version `2.0`. Please use the `cellrank.kernels` or `cellrank.estimators` interface instead.
  k = cr.tl.transition_matrix(

(0.8 * VelocityKernel + 0.2 * ConnectivityKernel)

Next, we compute the projection of the transition matrix onto an embedding, in this case UMAP. The projected transition matrix can be found in adata.obsm['T_fwd_umap'].

k.plot_projection(basis="umap")

The projected transition matrix can now be visualized using scvelo. Below we use scvelo.pl.velocity_embedding_stream() to visualize it as streamlines.

scv.pl.velocity_embedding_stream(adata, vkey="T_fwd", basis="umap")

Another way of visualizing the transition matrix is to simulate random walks on Markov chain and plot them in an embedding, see Plot random walks for more information.

Estimated memory usage: 534 MB