ValueError: 'x0' must only have one dimension. when do dyn.pd.least_action()

[dyinboisry4u]

Hi, Dynamo group I'm used dynamo to analysis my scRNA-Seq data, It runs well before the LAP analyses, but I stuck at dyn.pd.least_action when do LAP analysis here is my code (omit the plot):

# prepare data
loom = sc.read_h5ad("./data/loom.h5ad.gz")
cd8 = sc.read_h5ad("./data/cd8.h5ad.gz")
cd8_loom = loom[loom.obs_names.isin(cd8.obs_names)].copy()
cd8_loom.obsm['X_umap'] = cd8.obsm['X_umap']
cd8_loom.obs['res06_anno'] = cd8.obs['res06_anno']

# velo
preprocessor = dyn.pp.Preprocessor(cell_cycle_score_enable=True)
preprocessor.preprocess_adata(cd8_loom, recipe='monocle')
dyn.tl.dynamics(cd8_loom, assumption_mRNA="ss", model='stochastic', est_method='negbin', cores=128)
dyn.tl.neighbors(cd8_loom, n_pca_components=30, n_neighbors=30)
dyn.tl.cell_velocities(cd8_loom, method='pearson', other_kernels_dict={'transform': 'sqrt'})
dyn.tl.cell_wise_confidence(cd8_loom)
dyn.tl.confident_cell_velocities(cd8_loom, group='res06_anno', lineage_dict={'Tn': ['Tex']})

# vector field
dyn.vf.VectorField(cd8_loom, basis='umap', method='sparseVFC', M=1000, MaxIter=500, pot_curl_div=True, cores=128)
dyn.tl.cell_velocities(cd8_loom, basis='pca')
dyn.vf.VectorField(cd8_loom, basis='pca', cores=128)
dyn.vf.speed(cd8_loom, basis='pca')
dyn.vf.curl(cd8_loom, basis='umap')
dyn.vf.divergence(cd8_loom, basis='pca')
dyn.vf.acceleration(cd8_loom, basis='pca')
dyn.vf.curvature(cd8_loom, basis='pca')

# fate
progenitor = cd8_loom.obs_names[cd8_loom.obs.res06_anno.isin(['Tn'])]
dyn.pd.fate(cd8_loom, basis='umap', init_cells=progenitor, interpolation_num=100, direction='forward', inverse_transform=False, average=False, cores=128)

# fate animation
fig, ax = plt.subplots()
ax = dyn.pl.topography(cd8_loom, color='res06_anno', ax=ax, save_show_or_return='return')
dyn.mv.animate_fates(cd8_loom, color='res06_anno', basis='umap', n_steps=100, fig=fig, ax=ax,
                     save_show_or_return='save', save_kwargs={'filename': 'test.gif', 'writer': 'pillow'}, logspace=True, max_time=None)

# LAP
from dynamo.tools.utils import nearest_neighbors

Tn_cells = dyn.tl.select_cell(cd8_loom, "res06_anno", "Tn")
Tex_cells = dyn.tl.select_cell(cd8_loom, "res06_anno", "Tex")

fix_idx = [10, 17]
fixed_points = cd8_loom.uns['VecFld_umap']['Xss'][fix_idx]

Tn_cells_indices = nearest_neighbors(fixed_points[0], cd8_loom.obsm["X_umap"])
Tex_cells_indices = nearest_neighbors(fixed_points[1], cd8_loom.obsm["X_umap"])

dyn.tl.neighbors(cd8_loom, basis="umap", result_prefix="X_umap")
dyn.tl.cell_velocities(cd8_loom, method='cosine', add_transition_key='cosine_transition_matrix', enforce=True)

## get an error
dyn.dynamo_logger.main_silence()
transition_graph = {}
cell_type = ["Tn", "Tex"]
start_cell_indices = [
    Tn_cells_indices,
    Tex_cells_indices,
]
end_cell_indices = start_cell_indices
for i, start in enumerate(start_cell_indices):
    for j, end in enumerate(end_cell_indices):
        if start is not end:
            min_lap_t = True if i == 0 else False
            dyn.pd.least_action(
                cd8_loom,
                [cd8_loom.obs_names[start[0]][0]],
                [cd8_loom.obs_names[end[0]][0]],
                basis="umap",
                adj_key="X_umap_distances",
                min_lap_t=min_lap_t,
                EM_steps=2,
            )
            dyn.pl.least_action(cd8_loom, basis="umap")
            lap = dyn.pd.least_action(
                cd8_loom,
                [cd8_loom.obs_names[start[0]][0]],
                [cd8_loom.obs_names[end[0]][0]],
                basis="pca",
                adj_key="cosine_transition_matrix",
                min_lap_t=min_lap_t,
                EM_steps=2,
            )
            dyn.pl.kinetic_heatmap(
                cd8_loom,
                basis="pca",
                mode="lap",
                genes=cd8_loom.var_names[cd8_loom.var.use_for_transition],
                project_back_to_high_dim=True,
            )
            # The `GeneTrajectory` class can be used to output trajectories for any set of genes of interest
            gtraj = dyn.pd.GeneTrajectory(cd8_loom)
            gtraj.from_pca(lap.X, t=lap.t)
            gtraj.calc_msd()
            ranking = dyn.vf.rank_genes(cd8_loom, "traj_msd")

            print(start, "->", end)
            genes = ranking[:5]["all"].to_list()
            arr = gtraj.select_gene(genes)

            dyn.pl.multiplot(lambda k: [plt.plot(arr[k, :]), plt.title(genes[k])], np.arange(len(genes)))

            transition_graph[cell_type[i] + "->" + cell_type[j]] = {
                "lap": lap,
                "LAP_umap": cd8_loom.uns["LAP_umap"],
                "LAP_pca": cd8_loom.uns["LAP_pca"],
                "ranking": ranking,
                "gtraj": gtraj,
            }

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[23], line 19
     17 if start is not end:
     18     min_lap_t = True if i == 0 else False
---> 19     dyn.pd.least_action(
     20         cd8_loom,
     21         [cd8_loom.obs_names[start[0]][0]],
     22         [cd8_loom.obs_names[end[0]][0]],
     23         basis="umap",
     24         adj_key="X_umap_distances",
     25         min_lap_t=min_lap_t,
     26         EM_steps=2,
     27     )
     28     dyn.pl.least_action(cd8_loom, basis="umap")
     29     lap = dyn.pd.least_action(
     30         cd8_loom,
     31         [cd8_loom.obs_names[start[0]][0]],
   (...)
     36         EM_steps=2,
     37     )

File ~/software/install_pkg/miniconda3/envs/velo/lib/python3.11/site-packages/dynamo/prediction/least_action_path.py:555, in least_action(adata, init_cells, target_cells, init_states, target_states, paired, min_lap_t, elbow_method, num_t, basis, vf_key, vecfld, adj_key, n_points, init_paths, D, PCs, expr_func, add_key, **kwargs)
    550 path_ind += 1
    551 logger.info(
    552     "optimizing for least action path...",
    553     indent_level=2,
    554 )
--> 555 path_sol, dt_sol, action_opt = least_action_path(
    556     init_state, target_state, vf.func, vf.get_Jacobian(), n_points=n_points, init_path=init_path, D=D, **kwargs
    557 )
    559 n_points = len(path_sol)  # the actual #points due to arclength resampling
    561 if min_lap_t:

File ~/software/install_pkg/miniconda3/envs/velo/lib/python3.11/site-packages/dynamo/prediction/least_action_path.py:378, in least_action_path(start, end, vf_func, jac_func, n_points, init_path, D, dt_0, EM_steps)
    376 while EM_steps > 0:
    377     EM_steps -= 1
--> 378     path, dt, action_opt = lap_T(path, dt * len(path), vf_func, jac_func, D=D)
    380 return path, dt, action_opt

File ~/software/install_pkg/miniconda3/envs/velo/lib/python3.11/site-packages/dynamo/prediction/least_action_path.py:325, in lap_T(path_0, T, vf_func, jac_func, D)
    322 def jac(x):
    323     return action_grad_aux(x, vf_func, jac_func, dim, start=path_0[0], end=path_0[-1], D=D, dt=dt)
--> 325 sol_dict = minimize(fun, path_0[1:-1], jac=jac)
    326 path_sol = reshape_path(sol_dict["x"], dim, start=path_0[0], end=path_0[-1])
    328 # further optimization by varying dt

File ~/software/install_pkg/miniconda3/envs/velo/lib/python3.11/site-packages/scipy/optimize/_minimize.py:533, in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
    530 x0 = np.atleast_1d(np.asarray(x0))
    532 if x0.ndim != 1:
--> 533     raise ValueError("'x0' must only have one dimension.")
    535 if x0.dtype.kind in np.typecodes["AllInteger"]:
    536     x0 = np.asarray(x0, dtype=float)

ValueError: 'x0' must only have one dimension.

# session_info
-----
anndata             0.9.2
dynamo              1.3.2
ipykernel           6.25.1
ipywidgets          8.1.1
matplotlib          3.5.3
numpy               1.24.0
pandas              2.1.1
scanpy              1.9.3
scipy               1.11.2
session_info        1.0.0
-----

-----
IPython             8.14.0
jupyter_client      8.3.0
jupyter_core        5.3.1
-----
Python 3.11.4 | packaged by conda-forge | (main, Jun 10 2023, 18:08:17) [GCC 12.2.0]
Linux-3.10.0-1160.49.1.el7.x86_64-x86_64-with-glibc2.17
-----
Session information updated at 2023-09-27 20:55

What should I do, thanks!

[Sichao25]

Hi, thank you for reporting this issue. The update in Scipy causes this bug. Multiple dimension input is deprecated after SciPy 1.11.0. The easiest to solve this is to downgrade Scipy to an earlier version. Another way is to flatten path_0[1:-1] when calling minimize if you want to edit your local source code. We will fix this compatibility issue soon.

[dyinboisry4u]

Hi, thank you for reporting this issue. The update in Scipy causes this bug. Multiple dimension input is deprecated after SciPy 1.11.0. The easiest to solve this is to downgrade Scipy to an earlier version. Another way is to flatten path_0[1:-1] when calling minimize if you want to edit your local source code. We will fix this compatibility issue soon.

Thanks sichao, it works.

[github-actions[bot]]

This issue is stale because it has been open 90 days with no activity. Remove stale label or comment or this will be closed in 14 days