Least action path analysis with Dynamo

Driver analysis on scEU-seq organoid data using Dynamo’s LAP analysis

import itertools
import sys

import numpy as np
import pandas as pd

import dynamo as dyn
import scanpy as sc

sys.path.extend(["../../../", "."])
from paths import DATA_DIR  # isort: skip  # noqa: E402

|-----> setting visualization default mode in dynamo. Your customized matplotlib settings might be overritten.

/vol/storage/miniconda3/envs/dynamo-py39/lib/python3.9/site-packages/nxviz/__init__.py:18: UserWarning: 
nxviz has a new API! Version 0.7.3 onwards, the old class-based API is being
deprecated in favour of a new API focused on advancing a grammar of network
graphics. If your plotting code depends on the old API, please consider
pinning nxviz at version 0.7.3, as the new API will break your old code.

To check out the new API, please head over to the docs at
https://ericmjl.github.io/nxviz/ to learn more. We hope you enjoy using it!

(This deprecation message will go away in version 1.0.)

  warnings.warn(

General settings

sc.settings.verbosity = 3

Constants

TERMINAL_STATES = ["Enterocytes", "Enteroendocrine progenitors", "Goblet cells", "Paneth cells"]

Data loading

adata = sc.read(DATA_DIR / "sceu_organoid" / "processed" / "adata_dynamo-chase_and_pulse-2000features.h5ad")
adata

AnnData object with n_obs × n_vars = 3452 × 2000
    obs: 'experiment', 'time', 'cell_type', 'som_cluster_id', 'cell_type_merged', 'initial_size', 'n_counts', 'pass_basic_filter', 'ntr'
    var: 'ensum_id', 'gene_count_corr', 'means', 'dispersions', 'dispersions_norm', 'highly_variable', 'pass_basic_filter', 'use_for_pca', 'ntr', 'alpha', 'a', 'b', 'alpha_a', 'alpha_i', 'beta', 'p_half_life', 'gamma', 'half_life', 'cost', 'logLL', 'use_for_dynamics'
    uns: 'PCs', 'cell_type_colors', 'dynamics', 'neighbors', 'pca', 'pca_fit', 'pca_mean', 'pp', 'umap'
    obsm: 'X_pca', 'X_umap', 'X_umap_paper'
    varm: 'PCs'
    layers: 'M_n', 'M_nn', 'M_t', 'M_tn', 'M_tt', 'X_new', 'X_total', 'new', 'total', 'velocity_N', 'velocity_T'
    obsp: 'connectivities', 'distances', 'moments_con'

goblet_markers = (
    pd.read_csv(DATA_DIR / "sceu_organoid" / "processed" / "goblet_markers.csv")["Gene"].str.lower().tolist()
)

goblet_markers = adata.var_names[adata.var_names.str.lower().isin(goblet_markers)]

goblet_regulators = (
    pd.read_csv(DATA_DIR / "sceu_organoid" / "processed" / "goblet_regulators.csv")["Gene"].str.lower().tolist()
)

goblet_regulators = adata.var_names[adata.var_names.str.lower().isin(goblet_regulators)]

goblet_and_paneth_regulators = (
    pd.read_csv(DATA_DIR / "sceu_organoid" / "processed" / "goblet_and_paneth_regulators.csv")["Gene"]
    .str.lower()
    .tolist()
)

goblet_and_paneth_regulators = adata.var_names[adata.var_names.str.lower().isin(goblet_and_paneth_regulators)]

paneth_markers = (
    pd.read_csv(DATA_DIR / "sceu_organoid" / "processed" / "paneth_markers.csv")["Gene"].str.lower().tolist()
)

paneth_markers = adata.var_names[adata.var_names.str.lower().isin(paneth_markers)]

eec_markers = pd.read_csv(DATA_DIR / "sceu_organoid" / "processed" / "eec_markers.csv")["Gene"].str.lower().tolist()

eec_markers = adata.var_names[adata.var_names.str.lower().isin(eec_markers)]

eec_progenitor_markers = (
    pd.read_csv(DATA_DIR / "sceu_organoid" / "processed" / "eec_progenitor_markers.csv")["Gene"].str.lower().tolist()
)

eec_progenitor_markers = adata.var_names[adata.var_names.str.lower().isin(eec_progenitor_markers)]

enterocyte_markers = (
    pd.read_csv(DATA_DIR / "sceu_organoid" / "processed" / "enterocyte_markers.csv")["Gene"].str.lower().tolist()
)

enterocyte_markers = adata.var_names[adata.var_names.str.lower().isin(enterocyte_markers)]

enterocyte_progenitor_markers = (
    pd.read_csv(DATA_DIR / "sceu_organoid" / "processed" / "enterocyte_progenitor_markers.csv")["Gene"]
    .str.lower()
    .tolist()
)

enterocyte_progenitor_markers = adata.var_names[adata.var_names.str.lower().isin(enterocyte_progenitor_markers)]

Preprocessing

markers = {}
markers["Goblet cells"] = goblet_markers.union(goblet_regulators).union(goblet_and_paneth_regulators)
markers["Paneth cells"] = paneth_markers.union(goblet_and_paneth_regulators)
markers["Enteroendocrine progenitors"] = eec_markers.union(eec_progenitor_markers)
markers["Enterocytes"] = enterocyte_markers.union(enterocyte_progenitor_markers)

dyn.tl.reduceDimension(adata, layer="X_new", enforce=True)

|-----> retrive data for non-linear dimension reduction...
|-----> perform umap...
|-----> [dimension_reduction projection] in progress: 100.0000%
|-----> [dimension_reduction projection] finished [26.0535s]

Least action path analysis

dyn.tl.cell_velocities(adata, ekey="M_n", vkey="velocity_N", enforce=True)
dyn.vf.VectorField(adata, basis="umap")
dyn.vf.topography(adata)

|-----> retrive data for non-linear dimension reduction...
|-----? adata already have basis umap. dimension reduction umap will be skipped! 
set enforce=True to re-performing dimension reduction.
|-----> [dimension_reduction projection] in progress: 100.0000%
|-----> [dimension_reduction projection] finished [0.0022s]
|-----> 0 genes are removed because of nan velocity values.
|-----> [calculating transition matrix via pearson kernel with sqrt transform.] in progress: 100.0000%
|-----> [calculating transition matrix via pearson kernel with sqrt transform.] finished [2.4608s]
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%
|-----> [projecting velocity vector to low dimensional embedding] finished [0.8680s]
|-----> VectorField reconstruction begins...
|-----> Retrieve X and V based on basis: UMAP. 
        Vector field will be learned in the UMAP space.
|-----> Generating high dimensional grids and convert into a row matrix.
|-----> Learning vector field with method: sparsevfc.
|-----> [SparseVFC] begins...
|-----> Sampling control points based on data velocity magnitude...
|-----> [SparseVFC] in progress: 100.0000%
|-----> [SparseVFC] finished [22.5206s]
|-----------> current cosine correlation between input velocities and learned velocities is less than 0.6. Make a 1-th vector field reconstruction trial.
|-----> [SparseVFC] begins...
|-----> Sampling control points based on data velocity magnitude...
|-----> [SparseVFC] in progress: 100.0000%
|-----> [SparseVFC] finished [23.5246s]
|-----------> current cosine correlation between input velocities and learned velocities is less than 0.6. Make a 2-th vector field reconstruction trial.
|-----> [SparseVFC] begins...
|-----> Sampling control points based on data velocity magnitude...
|-----> [SparseVFC] in progress: 100.0000%
|-----> [SparseVFC] finished [23.9717s]
|-----------> current cosine correlation between input velocities and learned velocities is less than 0.6. Make a 3-th vector field reconstruction trial.
|-----> [SparseVFC] begins...
|-----> Sampling control points based on data velocity magnitude...
|-----> [SparseVFC] in progress: 100.0000%
|-----> [SparseVFC] finished [22.9951s]
|-----------> current cosine correlation between input velocities and learned velocities is less than 0.6. Make a 4-th vector field reconstruction trial.
|-----> [SparseVFC] begins...
|-----> Sampling control points based on data velocity magnitude...
|-----> [SparseVFC] in progress: 100.0000%
|-----> [SparseVFC] finished [22.1381s]
|-----------> current cosine correlation between input velocities and learned velocities is less than 0.6. Make a 5-th vector field reconstruction trial.
|-----? Cosine correlation between input velocities and learned velocities is less than 0.6 after 5 trials of vector field reconstruction.
|-----> <insert> velocity_umap_SparseVFC to obsm in AnnData Object.
|-----> <insert> X_umap_SparseVFC to obsm in AnnData Object.
|-----> <insert> VecFld_umap to uns in AnnData Object.
|-----> <insert> control_point_umap to obs in AnnData Object.
|-----> <insert> inlier_prob_umap to obs in AnnData Object.
|-----> <insert> obs_vf_angle_umap to obs in AnnData Object.
|-----> [VectorField] in progress: 100.0000%
|-----> [VectorField] finished [115.3775s]

/vol/storage/miniconda3/envs/dynamo-py39/lib/python3.9/site-packages/scipy/optimize/_minpack_py.py:178: RuntimeWarning: The iteration is not making good progress, as measured by the 
  improvement from the last ten iterations.
  warnings.warn(msg, RuntimeWarning)
/vol/storage/miniconda3/envs/dynamo-py39/lib/python3.9/site-packages/scipy/optimize/_minpack_py.py:178: RuntimeWarning: The iteration is not making good progress, as measured by the 
  improvement from the last five Jacobian evaluations.
  warnings.warn(msg, RuntimeWarning)

AnnData object with n_obs × n_vars = 3452 × 2000
    obs: 'experiment', 'time', 'cell_type', 'som_cluster_id', 'cell_type_merged', 'initial_size', 'n_counts', 'pass_basic_filter', 'ntr', 'control_point_umap', 'inlier_prob_umap', 'obs_vf_angle_umap'
    var: 'ensum_id', 'gene_count_corr', 'means', 'dispersions', 'dispersions_norm', 'highly_variable', 'pass_basic_filter', 'use_for_pca', 'ntr', 'alpha', 'a', 'b', 'alpha_a', 'alpha_i', 'beta', 'p_half_life', 'gamma', 'half_life', 'cost', 'logLL', 'use_for_dynamics', 'use_for_transition'
    uns: 'PCs', 'cell_type_colors', 'dynamics', 'neighbors', 'pca', 'pca_fit', 'pca_mean', 'pp', 'umap', 'explained_variance_ratio_', 'pca_valid_ind', 'X_new_neighbors', 'umap_fit', 'grid_velocity_umap', 'VecFld_umap'
    obsm: 'X_pca', 'X_umap', 'X_umap_paper', 'X_new_pca', 'X_new_umap', 'velocity_umap', 'velocity_umap_SparseVFC', 'X_umap_SparseVFC'
    varm: 'PCs'
    layers: 'M_n', 'M_nn', 'M_t', 'M_tn', 'M_tt', 'X_new', 'X_total', 'new', 'total', 'velocity_N', 'velocity_T'
    obsp: 'connectivities', 'distances', 'moments_con', 'pearson_transition_matrix'

adata.obsm["X_pca_orig"] = adata.obsm["X_pca"].copy()
adata.obsm["X_pca"] = adata.obsm["X_pca"][:, :30]

dyn.tl.cell_velocities(adata, ekey="M_n", vkey="velocity_N", basis="pca")
dyn.vf.VectorField(adata, basis="pca")

|-----> retrive data for non-linear dimension reduction...
|-----? adata already have basis pca. dimension reduction pca will be skipped! 
set enforce=True to re-performing dimension reduction.
|-----> [dimension_reduction projection] in progress: 100.0000%
|-----> [dimension_reduction projection] finished [0.0024s]
|-----> 0 genes are removed because of nan velocity values.
Using existing pearson_transition_matrix found in .obsp.
|-----> [projecting velocity vector to low dimensional embedding] in progress: 100.0000%
|-----> [projecting velocity vector to low dimensional embedding] finished [3.4971s]
|-----> VectorField reconstruction begins...
|-----> Retrieve X and V based on basis: PCA. 
        Vector field will be learned in the PCA space.
|-----> Learning vector field with method: sparsevfc.
|-----> [SparseVFC] begins...
|-----> Sampling control points based on data velocity magnitude...
|-----> [SparseVFC] in progress: 100.0000%
|-----> [SparseVFC] finished [3.6416s]
|-----> <insert> velocity_pca_SparseVFC to obsm in AnnData Object.
|-----> <insert> X_pca_SparseVFC to obsm in AnnData Object.
|-----> <insert> VecFld_pca to uns in AnnData Object.
|-----> <insert> control_point_pca to obs in AnnData Object.
|-----> <insert> inlier_prob_pca to obs in AnnData Object.
|-----> <insert> obs_vf_angle_pca to obs in AnnData Object.
|-----> [VectorField] in progress: 100.0000%
|-----> [VectorField] finished [3.8051s]

df = pd.DataFrame(adata.uns["VecFld_umap"]["Xss"], columns=["umap_1", "umap_2"])
df["fixed_point_type"] = adata.uns["VecFld_umap"]["ftype"]
df["fixed_point_type"].replace({-1: "stable", 0: "saddle", 1: "unstable"}, inplace=True)

neighbor_idx = dyn.tools.utils.nearest_neighbors(df[["umap_1", "umap_2"]], adata.obsm["X_new_umap"])
df["cell_type"] = [
    adata.obs.loc[adata.obs_names[neighbors], "cell_type"].mode().values[0] for neighbors in neighbor_idx
]

initial_obs = dyn.tools.utils.nearest_neighbors(
    df.loc[(df["fixed_point_type"] == "unstable") & (df["cell_type"] == "Stem cells"), ["umap_1", "umap_2"]].values,
    adata.obsm["X_new_umap"],
    k=30,
).flatten()

identified_terminal_states = list(
    set(TERMINAL_STATES).intersection(df.loc[df["fixed_point_type"] == "stable", "cell_type"].unique())
)

terminal_obs = {}
for terminal_state in identified_terminal_states:
    terminal_obs[terminal_state] = dyn.tools.utils.nearest_neighbors(
        df.loc[(df["fixed_point_type"] == "stable") & (df["cell_type"] == terminal_state), ["umap_1", "umap_2"]].values,
        adata.obsm["X_new_umap"],
        k=30,
    ).flatten()

dyn.tl.neighbors(adata, basis="umap", result_prefix="umap")

|-----> Start computing neighbor graph...
|-----------> X_data is None, fetching or recomputing...
|-----> fetching X data from layer:None, basis:umap
|-----> method arg is None, choosing methods automatically...
|-----------> method kd_tree selected
|-----> <insert> umap_connectivities to obsp in AnnData Object.
|-----> <insert> umap_distances to obsp in AnnData Object.
|-----> <insert> umap_neighbors to uns in AnnData Object.
|-----> <insert> umap_neighbors.indices to uns in AnnData Object.
|-----> <insert> umap_neighbors.params to uns in AnnData Object.

AnnData object with n_obs × n_vars = 3452 × 2000
    obs: 'experiment', 'time', 'cell_type', 'som_cluster_id', 'cell_type_merged', 'initial_size', 'n_counts', 'pass_basic_filter', 'ntr', 'control_point_umap', 'inlier_prob_umap', 'obs_vf_angle_umap', 'control_point_pca', 'inlier_prob_pca', 'obs_vf_angle_pca'
    var: 'ensum_id', 'gene_count_corr', 'means', 'dispersions', 'dispersions_norm', 'highly_variable', 'pass_basic_filter', 'use_for_pca', 'ntr', 'alpha', 'a', 'b', 'alpha_a', 'alpha_i', 'beta', 'p_half_life', 'gamma', 'half_life', 'cost', 'logLL', 'use_for_dynamics', 'use_for_transition'
    uns: 'PCs', 'cell_type_colors', 'dynamics', 'neighbors', 'pca', 'pca_fit', 'pca_mean', 'pp', 'umap', 'explained_variance_ratio_', 'pca_valid_ind', 'X_new_neighbors', 'umap_fit', 'grid_velocity_umap', 'VecFld_umap', 'grid_velocity_pca', 'VecFld_pca', 'umap_neighbors'
    obsm: 'X_pca', 'X_umap', 'X_umap_paper', 'X_new_pca', 'X_new_umap', 'velocity_umap', 'velocity_umap_SparseVFC', 'X_umap_SparseVFC', 'X_pca_orig', 'velocity_pca', 'velocity_pca_SparseVFC', 'X_pca_SparseVFC'
    varm: 'PCs'
    layers: 'M_n', 'M_nn', 'M_t', 'M_tn', 'M_tt', 'X_new', 'X_total', 'new', 'total', 'velocity_N', 'velocity_T'
    obsp: 'connectivities', 'distances', 'moments_con', 'pearson_transition_matrix', 'umap_distances', 'umap_connectivities'

Dynamo stable fixed points

gene_ranks = {}

for terminal_state in terminal_obs.keys():
    rankings = []

    lst = list(itertools.product(initial_obs, terminal_obs[terminal_state]))

    np.random.seed(0)
    state_tuples = np.random.choice(len(lst), size=10, replace=False).tolist()

    initial_obs_, terminal_obs_ = zip(*[lst[idx] for idx in state_tuples])

    laps = dyn.pd.least_action(
        adata,
        init_cells=list(initial_obs_),
        target_cells=list(terminal_obs_),
    )

    gtraj = dyn.pd.GeneTrajectory(adata)

    for lap_id, lap in enumerate(laps):
        gtraj.from_pca(lap.X, t=lap.t)
        gtraj.calc_msd()

        ranking = dyn.vf.rank_genes(adata, "traj_msd")
        ranking = ranking.reset_index().rename(columns={"index": f"Corr. rank - {terminal_state}", "all": "Gene"})
        ranking["Algorithm"] = "Dynamo"
        ranking["Run"] = lap_id
        rankings.append(ranking)

    gene_ranks[terminal_state] = pd.concat(rankings)
    gene_ranks[terminal_state] = gene_ranks[terminal_state].loc[
        gene_ranks[terminal_state]["Gene"].isin(markers[terminal_state])
    ]

|-----> searching for the least action path...
|-----> [iterating through 10 pairs] in progress: 10.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.043242
|-----> [iterating through 10 pairs] in progress: 20.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.048698
|-----> [iterating through 10 pairs] in progress: 30.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.036693
|-----> [iterating through 10 pairs] in progress: 40.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.054641
|-----> [iterating through 10 pairs] in progress: 50.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.036684
|-----> [iterating through 10 pairs] in progress: 60.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.040760
|-----> [iterating through 10 pairs] in progress: 70.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.052475
|-----> [iterating through 10 pairs] in progress: 80.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.055767
|-----> [iterating through 10 pairs] in progress: 90.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.081662
|-----> [iterating through 10 pairs] in progress: 100.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.043740
|-----> [iterating through 10 pairs] in progress: 100.0000%
|-----> [iterating through 10 pairs] finished [378.5261s]
|-----> [least action path] in progress: 100.0000%
|-----> [least action path] finished [378.5285s]
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.

/vol/storage/miniconda3/envs/dynamo-py39/lib/python3.9/site-packages/dynamo/prediction/trajectory.py:223: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  self.adata.var[save_key][self.genes_to_mask()] = msd

|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.

for terminal_state in gene_ranks.keys():
    gene_ranks[terminal_state].set_index("Gene", inplace=True)
    gene_ranks[terminal_state].index.name = None
    gene_ranks[terminal_state].to_csv(
        DATA_DIR
        / "sceu_organoid"
        / "results"
        / f"gene_ranks_{terminal_state}-chase_and_pulse-dynamo_terminal_states-dynamo.csv"
    )

CellRank terminal states

cr_terminal_states = pd.read_csv(
    DATA_DIR / "sceu_organoid" / "results" / "cr_terminal_states.csv", index_col=0
).reset_index(drop=True)["terminal_state"]
terminal_obs = {
    terminal_state: cr_terminal_states.index[cr_terminal_states == terminal_state].to_numpy()
    for terminal_state in cr_terminal_states.astype("category").cat.categories
}

gene_ranks = {}

for terminal_state in terminal_obs.keys():
    rankings = []

    lst = list(itertools.product(initial_obs, terminal_obs[terminal_state]))

    np.random.seed(0)
    state_tuples = np.random.choice(len(lst), size=10, replace=False).tolist()

    initial_obs_, terminal_obs_ = zip(*[lst[idx] for idx in state_tuples])

    laps = dyn.pd.least_action(
        adata,
        init_cells=list(initial_obs_),
        target_cells=list(terminal_obs_),
    )

    gtraj = dyn.pd.GeneTrajectory(adata)

    for lap_id, lap in enumerate(laps):
        gtraj.from_pca(lap.X, t=lap.t)
        gtraj.calc_msd()

        ranking = dyn.vf.rank_genes(adata, "traj_msd")
        ranking = ranking.reset_index().rename(columns={"index": f"Corr. rank - {terminal_state}", "all": "Gene"})
        ranking["Algorithm"] = "Dynamo"
        ranking["Run"] = lap_id
        rankings.append(ranking)

    gene_ranks[terminal_state] = pd.concat(rankings)
    gene_ranks[terminal_state] = gene_ranks[terminal_state].loc[
        gene_ranks[terminal_state]["Gene"].isin(markers[terminal_state])
    ]

|-----> searching for the least action path...
|-----> [iterating through 10 pairs] in progress: 10.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.056222
|-----> [iterating through 10 pairs] in progress: 20.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.044589
|-----> [iterating through 10 pairs] in progress: 30.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.044448
|-----> [iterating through 10 pairs] in progress: 40.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.058033
|-----> [iterating through 10 pairs] in progress: 50.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.046548
|-----> [iterating through 10 pairs] in progress: 60.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.051713
|-----> [iterating through 10 pairs] in progress: 70.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.087245
|-----> [iterating through 10 pairs] in progress: 80.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.062437
|-----> [iterating through 10 pairs] in progress: 90.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.042019
|-----> [iterating through 10 pairs] in progress: 100.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.053537
|-----> [iterating through 10 pairs] in progress: 100.0000%
|-----> [iterating through 10 pairs] finished [356.4280s]
|-----> [least action path] in progress: 100.0000%
|-----> [least action path] finished [356.4302s]
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.

/vol/storage/miniconda3/envs/dynamo-py39/lib/python3.9/site-packages/dynamo/prediction/trajectory.py:223: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  self.adata.var[save_key][self.genes_to_mask()] = msd

|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> searching for the least action path...
|-----> [iterating through 10 pairs] in progress: 10.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.039659
|-----> [iterating through 10 pairs] in progress: 20.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.029670
|-----> [iterating through 10 pairs] in progress: 30.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.031246
|-----> [iterating through 10 pairs] in progress: 40.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.024404
|-----> [iterating through 10 pairs] in progress: 50.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.036357
|-----> [iterating through 10 pairs] in progress: 60.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.032059
|-----> [iterating through 10 pairs] in progress: 70.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.026341
|-----> [iterating through 10 pairs] in progress: 80.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.026243
|-----> [iterating through 10 pairs] in progress: 90.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.032776
|-----> [iterating through 10 pairs] in progress: 100.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.022523
|-----> [iterating through 10 pairs] in progress: 100.0000%
|-----> [iterating through 10 pairs] finished [265.2149s]
|-----> [least action path] in progress: 100.0000%
|-----> [least action path] finished [265.2174s]
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.

/vol/storage/miniconda3/envs/dynamo-py39/lib/python3.9/site-packages/dynamo/prediction/trajectory.py:223: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  self.adata.var[save_key][self.genes_to_mask()] = msd

|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> searching for the least action path...
|-----> [iterating through 10 pairs] in progress: 10.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.040390
|-----> [iterating through 10 pairs] in progress: 20.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.036268
|-----> [iterating through 10 pairs] in progress: 30.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.039184
|-----> [iterating through 10 pairs] in progress: 40.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.026063
|-----> [iterating through 10 pairs] in progress: 50.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.028168
|-----> [iterating through 10 pairs] in progress: 60.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.014053
|-----> [iterating through 10 pairs] in progress: 70.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.042330
|-----> [iterating through 10 pairs] in progress: 80.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.025680
|-----> [iterating through 10 pairs] in progress: 90.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.020968
|-----> [iterating through 10 pairs] in progress: 100.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.026264
|-----> [iterating through 10 pairs] in progress: 100.0000%
|-----> [iterating through 10 pairs] finished [236.9429s]
|-----> [least action path] in progress: 100.0000%
|-----> [least action path] finished [236.9447s]
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.

/vol/storage/miniconda3/envs/dynamo-py39/lib/python3.9/site-packages/dynamo/prediction/trajectory.py:223: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  self.adata.var[save_key][self.genes_to_mask()] = msd

|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> searching for the least action path...
|-----> [iterating through 10 pairs] in progress: 10.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.024865
|-----> [iterating through 10 pairs] in progress: 20.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.024610
|-----> [iterating through 10 pairs] in progress: 30.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.019792
|-----> [iterating through 10 pairs] in progress: 40.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.021137
|-----> [iterating through 10 pairs] in progress: 50.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.022025
|-----> [iterating through 10 pairs] in progress: 60.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.023605
|-----> [iterating through 10 pairs] in progress: 70.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.022567
|-----> [iterating through 10 pairs] in progress: 80.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.042867
|-----> [iterating through 10 pairs] in progress: 90.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.024776
|-----> [iterating through 10 pairs] in progress: 100.0000%|-----------> initializing path with the shortest path in the graph built from the velocity transition matrix...
|-----------> optimizing for least action path...
|-----> optimal action: 0.013294
|-----> [iterating through 10 pairs] in progress: 100.0000%
|-----> [iterating through 10 pairs] finished [281.3232s]
|-----> [least action path] in progress: 100.0000%
|-----> [least action path] finished [281.3252s]
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.

/vol/storage/miniconda3/envs/dynamo-py39/lib/python3.9/site-packages/dynamo/prediction/trajectory.py:223: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  self.adata.var[save_key][self.genes_to_mask()] = msd

|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.
|-----> <insert> traj_msd to var in AnnData Object.

for terminal_state in gene_ranks.keys():
    gene_ranks[terminal_state].set_index("Gene", inplace=True)
    gene_ranks[terminal_state].index.name = None
    gene_ranks[terminal_state].to_csv(
        DATA_DIR
        / "sceu_organoid"
        / "results"
        / f"gene_ranks_{terminal_state}-chase_and_pulse-cr_terminal_states-dynamo.csv"
    )