Is there a way to run the animate fate transition code block for an in-silico perturbed vector field?

[kousaa]

Hi Dynamo team,

thank you so much for developing this great tool. I have been playing around with the in-silico perturbation function and I came to wonder if it would be possible to create a gif animation of the potential trajectory on the perturbed space.

I tried to reconstruct the vector field based on the perturbed UMAP dyn.vf.VectorField(adata, basis='umap_perturbation', M=1000, MaxIter=170, pot_curl_div=True)

and then run the animation function

dyn.pd.fate(adata, basis='umap_perturbation', init_cells=progenitor, interpolation_num=100, direction='forward', inverse_transform=False, average=False, cores=7)

but I get the following error

calc_fft_k run failed...
calc_fft_k run failed...
calc_fft_k run failed...
calc_fft_k run failed...
calc_fft_k run failed...

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[166], line 1
----> 1 dyn.pd.fate(adata, basis='umap_perturbation', init_cells=progenitor, interpolation_num=100,  direction='forward',
      2    inverse_transform=False, average=False, cores=7)

File ~/.local/lib/python3.8/site-packages/dynamo/prediction/fate.py:136, in fate(adata, init_cells, init_states, basis, layer, dims, genes, t_end, direction, interpolation_num, average, sampling, VecFld_true, inverse_transform, Qkey, scale, cores, **kwargs)
    131     init_states = init_states[:, dims]
    133 vf = (
    134     (lambda x: scale * vector_field_function(x=x, vf_dict=VecFld, dim=dims)) if VecFld_true is None else VecFld_true
    135 )
--> 136 t, prediction = _fate(
    137     vf,
    138     init_states,
    139     t_end=t_end,
    140     direction=direction,
    141     interpolation_num=interpolation_num,
    142     average=True if average in ["origin", "trajectory", True] else False,
    143     sampling=sampling,
    144     cores=cores,
    145     **kwargs,
    146 )
    148 exprs = None
    149 if basis == "pca" and inverse_transform:

File ~/.local/lib/python3.8/site-packages/dynamo/prediction/fate.py:260, in _fate(VecFld, init_states, t_end, step_size, direction, interpolation_num, average, sampling, cores)
    258 else:
    259     pool = ThreadPool(cores)
--> 260     res = pool.starmap(
    261         integrate_vf_ivp,
    262         zip(
    263             init_states,
    264             itertools.repeat(t_linspace),
    265             itertools.repeat(direction),
    266             itertools.repeat(VecFld),
    267             itertools.repeat(()),
    268             itertools.repeat(interpolation_num),
    269             itertools.repeat(average),
    270             itertools.repeat(sampling),
    271             itertools.repeat(False),
    272             itertools.repeat(True),
    273         ),
    274     )  # disable tqdm when using multiple cores.
    275     pool.close()
    276     pool.join()

File ~/opt/anaconda3/envs/Dynamo/lib/python3.8/multiprocessing/pool.py:372, in Pool.starmap(self, func, iterable, chunksize)
    366 def starmap(self, func, iterable, chunksize=None):
    367     '''
    368     Like `map()` method but the elements of the `iterable` are expected to
    369     be iterables as well and will be unpacked as arguments. Hence
    370     `func` and (a, b) becomes func(a, b).
    371     '''
--> 372     return self._map_async(func, iterable, starmapstar, chunksize).get()

File ~/opt/anaconda3/envs/Dynamo/lib/python3.8/multiprocessing/pool.py:771, in ApplyResult.get(self, timeout)
    769     return self._value
    770 else:
--> 771     raise self._value

File ~/opt/anaconda3/envs/Dynamo/lib/python3.8/multiprocessing/pool.py:125, in worker(inqueue, outqueue, initializer, initargs, maxtasks, wrap_exception)
    123 job, i, func, args, kwds = task
    124 try:
--> 125     result = (True, func(*args, **kwds))
    126 except Exception as e:
    127     if wrap_exception and func is not _helper_reraises_exception:

File ~/opt/anaconda3/envs/Dynamo/lib/python3.8/multiprocessing/pool.py:51, in starmapstar(args)
     50 def starmapstar(args):
---> 51     return list(itertools.starmap(args[0], args[1]))

File ~/.local/lib/python3.8/site-packages/dynamo/prediction/utils.py:178, in integrate_vf_ivp(init_states, t, integration_direction, f, args, interpolation_num, average, sampling, verbose, disable)
    174     if integration_direction == "both":
    175         neg_t_len = sum(np.array(t_[i]) < 0)
    177     odeint_cur_Y = (
--> 178         SOL[i](t_[i])
    179         if integration_direction != "both"
    180         else np.hstack(
    181             (
    182                 SOL[i][0](t_[i][:neg_t_len]),
    183                 SOL[i][1](t_[i][neg_t_len:]),
    184             )
    185         )
    186     )
    187     Y_[i] = odeint_cur_Y
    189 Y, t = Y_, t_

File ~/opt/anaconda3/envs/Dynamo/lib/python3.8/site-packages/scipy/integrate/_ivp/common.py:236, in OdeSolution.__call__(self, t)
    233     ys.append(y)
    234     group_start = group_end
--> 236 ys = np.hstack(ys)
    237 ys = ys[:, reverse]
    239 return ys

File <__array_function__ internals>:200, in hstack(*args, **kwargs)

File ~/opt/anaconda3/envs/Dynamo/lib/python3.8/site-packages/numpy/core/shape_base.py:370, in hstack(tup, dtype, casting)
    368     return _nx.concatenate(arrs, 0, dtype=dtype, casting=casting)
    369 else:
--> 370     return _nx.concatenate(arrs, 1, dtype=dtype, casting=casting)

File <__array_function__ internals>:200, in concatenate(*args, **kwargs)

ValueError: need at least one array to concatenate

If you could help, that would be brilliant!

All the best, Anastasia

package 	umap-learn 	typing-extensions 	tqdm 	statsmodels 	setuptools 	session-info 	seaborn 	scipy 	scikit-learn 	pynndescent 	pre-commit 	pandas 	openpyxl 	numpy 	numdifftools 	numba 	networkx 	matplotlib 	loompy 	igraph 	get-version 	colorcet 	anndata 	dynamo-release
version 	0.5.3 	4.6.3 	4.65.0 	0.14.0 	67.8.0 	1.0.0 	0.12.2 	1.10.1 	1.3.0 	0.5.10 	3.3.3 	2.0.3 	3.1.2 	1.24.4 	0.9.41 	0.57.1 	3.1 	3.7.1 	3.0.7 	0.10.5 	3.5.4 	3.0.1 	0.9.1 	1.3.3

[Sichao25]

Hi, thanks for using Dynamo. We appreciate your feedback and contribution.

Just from the error message, it looks like Scipy fails to solve the initial value problem correctly. Your code looks good to me, so I suspect that the issue may be related to the specifics of the input data. Are you working with one of the sample datasets? Could you share your code related to data processing? If you can provide more information, it will help me better figure out the problem.

[kousaa]

Hi, thanks for trying to assist with this.

I can share the code for sure, but just before that I should say that the fate animation code block works just fine for my normal umap, it is only when I try to apply it on the perturbed one that it fails with that error. Do you think it would be possible to try to run the fate animation with a perturbed version of the example data you provide? If that works for you and you share the functional code, I could try to debug it on my end, figure out what goes wrong and share it here for anyone who may have the same issue. Thank you!

[Sichao25]

I see. At first, I successfully run the code on perturbed sample data. After spending some time fine-tuning the parameters, I reproduced the error you mentioned. It seems that the error is caused by a super short trajectory predicted in fate(). While one possibility is a bug in the implementation, it's more likely a data-specific special case that I can address by excluding that particular trajectory. I probably need a little more time to investigate this before I make any modifications to the code. Thanks again for raising this interesting issue.

You may be able to make a workaround by tuning parameters (like M, MaxIter, interpolation_num...), as I did manage to achieve success with the sample data. I will keep you updated with the final solution shortly.

[kousaa]

Thanks so much for your help with this! That makes a lot of sense; I will go ahead and tweak the tuning parameters so it works on my end too. :-)

A quick last question about this would be if there is a way to quantify the GIF animation in a static image. We are submitting a rebuttal for our paper and even though we will be including the GIFs as a movie, it would be great to have a visualization to show on the main figure of how many cells starting from cell A end up in the different subsets.

Would you have any suggestions?

Again, thanks so much for all your help and also developing this great tool.

[Sichao25]

When you mention "quantifying the GIF in a static image," are you referring to generating an image showing all the trajectories present in the animation? The GIF itself is actually an animated visualization of fate dynamics. We do have standard matplotlib plot functions for fate like .pl.fate(), which does something like:

ax = scatters(adata, basis=basis, color=color, save_show_or_return="return", ax=ax, **kwargs)
fate_key = "fate" if basis is None else "fate_" + basis
lap_dict = adata.uns[fate_key]
for i, j in zip(lap_dict["prediction"], lap_dict["t"]):
    ax.scatter(*i.T[:, [x, y]].T, c=map2color(j))
    ax.plot(*i.T[:, [x, y]].T, c="k")

You can modify the lap_dict here to select the trajectories you want to visualize. Also, the least action path is a good alternative for this. Let me know if this helps you.

[kousaa]

Thanks for the quick response! Let me try to explain myself better. Ideally, what I need is a visualization to show starting from population A what is the relative frequency of cells that are predicted to reach each of the other subsets.

There is a nice example visualization on the graphical abstract of this paper, where t for us would be probably be ddodge potential.

https://www.cell.com/cell/fulltext/S0092-8674(21)00439-6?_returnURL=https%3A%2F%2Flinkinghub.elsevier.com%2Fretrieve%2Fpii%2FS0092867421004396%3Fshowall%3Dtrue

But even if I could extract somehow the proportions of cell type A that become cell types B, C and E would suffice!

Thanks again!

[kousaa]

Hello!

I have being toying around in an effort to generate the above plot. I wanted to share my progress with you and test if this looks legit on your end.

I have used some of the functionality provided in CellRank (https://cellrank.readthedocs.io/en/latest/notebooks/tutorials/general/100_getting_started.html) to create a kernel and compute a transition matrix based on Dynamo's ddhodge potential.

from cellrank.kernels import PseudotimeKernel

dd = PseudotimeKernel(adata, time_key="umap_ddhodge_potential")
dd.compute_transition_matrix()

dd.plot_random_walks(
    seed=0,
    n_sims=150,
    start_ixs={"cell_type_subset": "subset A"},
    basis="umap",
    legend_loc="right",
    dpi=150,
)

Next, in order to use CellRanks's plot_single_flow function to create the same plot across ddhodge potential, I first had to bin the ddhodge potential as I show below, and then the plotting function worked just fine.

bins = [0,1,2,3,4,5,6,7]
adata.obs["umap_ddhodge_potential_binned"] = np.searchsorted(bins, adata.obs["umap_ddhodge_potential"].tolist())

ax = dd.plot_single_flow(
    cluster_key="cell_type_subset",
    time_key="umap_ddhodge_potential_binned",
    cluster="subset A",
    xticks_step_size=1,
    show=False,
        
    clusters=[ "subset B", "subset C", "subset D"],
)

_ = ax.set_xticklabels(ax.get_xticklabels(), rotation=90)

Now, this all works well, but I wanted to get a second opinion if this looks legit. I was also trying to use the pearson or cosine transition matrix from Dynamo instead calculating a transition metrix using the random walk algorithm, but when I tried to include it in the kernel I got the following error.

dd.transition_matrix = adata.obsp['pearson_transition_matrix']

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[129], line 1
----> 1 dd.transition_matrix = adata.obsp['pearson_transition_matrix']

File ~/opt/anaconda3/envs/Dynamo_CellRank/lib/python3.8/site-packages/cellrank/kernels/_base_kernel.py:595, in Kernel.transition_matrix(self, matrix)
    593 @transition_matrix.setter
    594 def transition_matrix(self, matrix: Any) -> None:
--> 595     KernelExpression.transition_matrix.fset(self, matrix)

File ~/opt/anaconda3/envs/Dynamo_CellRank/lib/python3.8/site-packages/cellrank/kernels/_base_kernel.py:429, in KernelExpression.transition_matrix(self, matrix)
    427 if force_normalize:
    428     if np.any((matrix.data if sp.issparse(matrix) else matrix) < 0):
--> 429         raise ValueError("Unable to normalize matrix with negative values.")
    430     matrix = _normalize(matrix)
    431     if should_norm(matrix):  # some rows are all 0s/contain invalid values

ValueError: Unable to normalize matrix with negative values.

Any advice or comments based on your expertise would be really appreciated.

Cheers, Anastasia

[Sichao25]

Hi. Sorry for the delayed response. I haven't checked the example visualization you provided yet so I will follow up on that later. About this specific problem, our transition matrix may contain negative values, which could be causing the error you encountered. I think what you can try is to set the use_neg_vals=False when you calculate the Pearson_transition_matrix with cell_velocities().

By the way, I think the original issue on the fate() will be fixed by a recent update. I can let you know when it is merged.

[kousaa]

Hello. thanks so much for this and sorry for the delayed response. Setting the use_neg_vals=False resolved my issue! Would you have an estimate for when the fate() fix update will be released?

Thanks again! Anastasia

[Sichao25]

It has been merged yesterday. You can try with the latest code now.

[kousaa]

Tested it and it works! Thank you so much for the help with this.

You can go ahead and close this issue as resolved, and I will create a new thread for anything else. :-)

Thanks again, Anastasia

[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