Pancreatic endocrinogenesis
[1]:
import math
import scanpy as sc
import scvelo as scv
import matplotlib.pyplot as plt
import sys
sys.path.append("..")
from tivelo.main import tivelo
from baseline import run_baseline
(Running UniTVelo 0.2.5.2)
2025-04-05 05:03:13
2025-04-05 13:03:14.006909: E tensorflow/stream_executor/cuda/cuda_blas.cc:2981] Unable to register cuBLAS factory: Attempting to register factory for plugin cuBLAS when one has already been registered
Run TIVelo
Load the dataset. Set the data name with key for cluster, key for embedding and cluster edges (for comparison).
[2]:
data_name = "pancreas"
data_path = "/lustre/project/Stat/s1155184322/datasets/velocity/endocrinogenesis_day15_processed.h5ad"
adata = sc.read(data_path)
group_key = "clusters"
emb_key = "X_umap"
cluster_edges = [("Pre-endocrine", "Alpha"), ("Pre-endocrine", "Beta"), ("Pre-endocrine", "Delta"),
("Pre-endocrine", "Epsilon")]
Set the model parameters.
[6]:
save_folder = "results"
show_fig = True
filter_genes = True
save_coeff = True
constrain = True
loss_fun = "mse"
only_s = False
alpha_1 = 1
alpha_2 = 0.1
batch_size = 1024
n_epochs = 100
tree_gene = None
show_DTI = False
adjust_DTI = False
velocity_key = "velocity"
measure_performance = True
Run the model by function tivelo.
[7]:
adata_ = tivelo(adata, group_key, emb_key, data_name=data_name, save_folder=save_folder, njobs=-1, tree_gene=tree_gene,
show_fig=show_fig, filter_genes=filter_genes, constrain=constrain, loss_fun=loss_fun, only_s=only_s,
alpha_1=alpha_1, alpha_2=alpha_2, batch_size=batch_size, n_epochs=n_epochs, velocity_key="velocity",
adjust_DTI=adjust_DTI, show_DTI=show_DTI, cluster_edges=cluster_edges,
measure_performance=measure_performance)
computing velocities
finished (0:00:00) --> added
'velocity', velocity vectors for each individual cell (adata.layers)
computing velocity graph (using 20/20 cores)
finished (0:00:13) --> added
'velocity_graph', sparse matrix with cosine correlations (adata.uns)
computing terminal states
WARNING: Uncertain or fuzzy root cell identification. Please verify.
identified 2 regions 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)
'path_dict' added to adata.uns
main path: ['Ductal', 'Ngn3 low EP', 'Ngn3 high EP', 'Pre-endocrine', 'Epsilon', 'Alpha', 'Beta']
[Parallel(n_jobs=-1)]: Using backend LokyBackend with 20 concurrent workers.
[Parallel(n_jobs=-1)]: Done 1 tasks | elapsed: 1.8s
[Parallel(n_jobs=-1)]: Done 10 tasks | elapsed: 2.1s
[Parallel(n_jobs=-1)]: Done 21 tasks | elapsed: 3.4s
[Parallel(n_jobs=-1)]: Done 32 tasks | elapsed: 3.9s
[Parallel(n_jobs=-1)]: Done 45 tasks | elapsed: 5.5s
[Parallel(n_jobs=-1)]: Done 58 tasks | elapsed: 6.0s
[Parallel(n_jobs=-1)]: Done 73 tasks | elapsed: 7.6s
[Parallel(n_jobs=-1)]: Done 88 tasks | elapsed: 9.2s
[Parallel(n_jobs=-1)]: Done 105 tasks | elapsed: 10.8s
[Parallel(n_jobs=-1)]: Done 122 tasks | elapsed: 12.4s
[Parallel(n_jobs=-1)]: Done 141 tasks | elapsed: 14.1s
[Parallel(n_jobs=-1)]: Done 160 tasks | elapsed: 15.2s
[Parallel(n_jobs=-1)]: Done 181 tasks | elapsed: 17.5s
[Parallel(n_jobs=-1)]: Done 202 tasks | elapsed: 19.2s
[Parallel(n_jobs=-1)]: Done 225 tasks | elapsed: 21.3s
[Parallel(n_jobs=-1)]: Done 248 tasks | elapsed: 23.5s
[Parallel(n_jobs=-1)]: Done 273 tasks | elapsed: 25.8s
[Parallel(n_jobs=-1)]: Done 298 tasks | elapsed: 28.3s
[Parallel(n_jobs=-1)]: Done 325 tasks | elapsed: 30.6s
[Parallel(n_jobs=-1)]: Done 352 tasks | elapsed: 32.9s
[Parallel(n_jobs=-1)]: Done 381 tasks | elapsed: 35.7s
[Parallel(n_jobs=-1)]: Done 410 tasks | elapsed: 38.2s
[Parallel(n_jobs=-1)]: Done 441 tasks | elapsed: 41.2s
[Parallel(n_jobs=-1)]: Done 472 tasks | elapsed: 43.8s
[Parallel(n_jobs=-1)]: Done 505 tasks | elapsed: 46.9s
[Parallel(n_jobs=-1)]: Done 538 tasks | elapsed: 50.0s
[Parallel(n_jobs=-1)]: Done 573 tasks | elapsed: 53.1s
[Parallel(n_jobs=-1)]: Done 608 tasks | elapsed: 56.1s
[Parallel(n_jobs=-1)]: Done 645 tasks | elapsed: 59.7s
[Parallel(n_jobs=-1)]: Done 682 tasks | elapsed: 1.0min
[Parallel(n_jobs=-1)]: Done 721 tasks | elapsed: 1.1min
[Parallel(n_jobs=-1)]: Done 760 tasks | elapsed: 1.2min
[Parallel(n_jobs=-1)]: Done 801 tasks | elapsed: 1.2min
[Parallel(n_jobs=-1)]: Done 842 tasks | elapsed: 1.3min
[Parallel(n_jobs=-1)]: Done 885 tasks | elapsed: 1.3min
[Parallel(n_jobs=-1)]: Done 928 tasks | elapsed: 1.4min
[Parallel(n_jobs=-1)]: Done 973 tasks | elapsed: 1.5min
[Parallel(n_jobs=-1)]: Done 1018 tasks | elapsed: 1.5min
[Parallel(n_jobs=-1)]: Done 1070 out of 1070 | elapsed: 1.6min finished
mean: 16.072
median: 10.323
lower quantile: -7.750
upper quantile: 35.858
minimum: -202.064
maximum: 312.487
No. of positive scores: 692
'path_dict' added to adata.uns
'child_dict' added to adata.uns
'level_dict' added to adata.uns
'threshold_list' added to adata.uns
'd_nn' added to adata.obsp
Model training: 100%|██████████| 100/100 [00:39<00:00, 2.52it/s, cos_s=0.623, cos_u=0.842, mse_s=0.008, mse_u=0.008]
computing velocity graph (using 20/20 cores)
finished (0:00:03) --> 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)
TIVelo:
CBDir: 0.5257 ICVCoh: 0.6353
CBDir2: 0.4850 ICVCoh2: 0.8043
TransProbs: 0.4745 VeloCoh: 0.1696
Run scVelo
Stochastic mode.
[8]:
adata_scvelo = run_baseline(adata, "scvelo", data_name, group_key, emb_key, cluster_edges, show_fig=True,
measure_performance=True)
computing velocity embedding
finished (0:00:00) --> added
'velocity_umap', embedded velocity vectors (adata.obsm)
scvelo:
CBDir: 0.5071 ICVCoh: 0.8699
CBDir2: 0.0413 ICVCoh2: 0.8081
TransProbs: 0.1732 VeloCoh: 0.0523
Dynamical mode.
[9]:
adata_scvelo2 = run_baseline(adata, "scvelo2", data_name, group_key, emb_key, cluster_edges, show_fig=True,
measure_performance=True)
scvelo2:
CBDir: 0.4430 ICVCoh: 0.8721
CBDir2: 0.1783 ICVCoh2: 0.8301
TransProbs: 0.2338 VeloCoh: 0.1971
Run veloVI.
[10]:
adata_velovi = run_baseline(adata, "velovi", data_name, group_key, emb_key, cluster_edges, show_fig=True,
measure_performance=True)
velovi:
CBDir: 0.5323 ICVCoh: 0.8940
CBDir2: 0.1042 ICVCoh2: 0.8271
TransProbs: 0.1065 VeloCoh: 0.1731
Run UniTVelo.
[3]:
adata_unitvelo = run_baseline(adata, "unitvelo", data_name, group_key, emb_key, cluster_edges, show_fig=True,
measure_performance=True, unitvelo_mode="2")
unitvelo:
CBDir: 0.4772 ICVCoh: 0.6340
CBDir2: 0.2001 ICVCoh2: 0.6824
TransProbs: 0.1953 VeloCoh: 0.1424
Run cellDancer
We don’t directly run cellDancer here since the environment required by cellDancer may conflict that of TIVelo. We recommend running cellDancer in an independent environment.
[12]:
adata_celldancer = run_baseline(adata, "celldancer", data_name, group_key, emb_key, cluster_edges, show_fig=True,
measure_performance=True)
computing velocity graph (using 20/20 cores)
finished (0:00:11) --> added
'velocity_S_graph', sparse matrix with cosine correlations (adata.uns)
computing velocity embedding
finished (0:00:00) --> added
'velocity_S_umap', embedded velocity vectors (adata.obsm)
celldancer:
CBDir: 0.4454 ICVCoh: 0.9548
CBDir2: 0.0583 ICVCoh2: 0.6712
TransProbs: 0.0479 VeloCoh: -0.2845
Run DeepVelo
We don’t directly run DeepVelo here since the environment required by DeepVelo may conflict that of TIVelo. We recommend running DeepVelo in an independent environment.
[3]:
adata_deepvelo = sc.read("/users/s1155184322/projects/tutorial/DeepVelo/results/pancreas_deepvelo.h5ad")
velocity_key = "velocity"
[4]:
ax = scv.pl.velocity_embedding_stream(adata_deepvelo, vkey=velocity_key, color=group_key, title="", show=False)
plt.tight_layout()
[5]:
from tivelo.utils.metrics import inner_cluster_coh, cross_boundary_correctness, cross_boundary_scvelo_probs, \
cross_boundary_correctness2, inner_cluster_coh2, velo_coh
_, cbdir = cross_boundary_correctness(adata_deepvelo, cluster_key=group_key, velocity_key=velocity_key,
cluster_edges=cluster_edges, x_emb=emb_key)
_, cbdir2 = cross_boundary_correctness2(adata_deepvelo, cluster_key=group_key, velocity_key=velocity_key,
cluster_edges=cluster_edges)
_, trans_probs = cross_boundary_scvelo_probs(adata_deepvelo, cluster_key=group_key, cluster_edges=cluster_edges,
trans_g_key="{}_graph".format(velocity_key))
_, icvcoh = inner_cluster_coh(adata_deepvelo, cluster_key=group_key, velocity_key=velocity_key)
_, icvcoh2 = inner_cluster_coh2(adata_deepvelo, cluster_key=group_key, velocity_key=velocity_key, x_emb=emb_key)
velocoh = velo_coh(adata_deepvelo, velocity_key=velocity_key, trans_g_key="{}_graph".format(velocity_key))
print("DeepVelo:\n", "CBDir:", "%.4f" % cbdir, "ICVCoh:", "%.4f" % icvcoh, "\n",
"CBDir2:", "%.4f" % cbdir2, "ICVCoh2:", "%.4f" % icvcoh2, "\n",
"TransProbs:", "%.4f" % trans_probs, "VeloCoh:", "%.4f" % velocoh)
DeepVelo:
CBDir: 0.4316 ICVCoh: 0.8947
CBDir2: 0.0074 ICVCoh2: 0.7977
TransProbs: 0.1097 VeloCoh: 0.1388