Plot directed PAGA

This example shows how to compute and plot a directed version of the popular PAGA algorithm [Wolf19].

In classical PAGA plot, nodes correspond to clusters and edge thickness denotes transcriptomic similarity. We introduce a new directed version of PAGA where directed edges reflect local velocity flow. We add the possibility to include prior information in the form of a pseudotemporal ordering or initial/terminal state annotation to restrict the possible edge set. We further replace nodes by pie charts that show average CellRank fate probabilities.

import cellrank as cr
import scvelo as scv

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 compute initial and terminal state probabilities as well as the absorption probabilities towards the identified terminal states.

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

cr.tl.initial_states(adata, cluster_key="clusters", n_states=1, softmax_scale=4)

Out:

INFO: Using pre-computed schur decomposition
HBox(children=(HTML(value=''), FloatProgress(value=0.0, max=2531.0), HTML(value='')))

We can use scvelo.tl.recover_latent_time() to compute gene-shared latent time leveraging the initial and terminal states computed above. This will be used as a time prior when computing the directed PAGA graph.

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

Out:

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

Afterwards, we compute the directed PAGA using scvelo.tl.paga() by again specifying the initial and terminal states and the time prior computed above.

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

Out:

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)

Lastly, we can plot the results using cellrank.pl.cluster_fates(). For more option, see Plot aggregated cellular fates.

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

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

Estimated memory usage: 679 MB

Gallery generated by Sphinx-Gallery