Troubleshooting issues with fitted model using scv.tl.recover_dynamics

[YitengDang]

Basically, I'm unable to fit the splicing dynamics to my data using scv.tl.recover_dynamics.

The fitting of the top likelihood genes seem to be very far off:

Some marker genes are fit slightly better, but also not as one would expect. Meanwhile, the stochastic model does seem to get the steady state of the marker genes right.

It seems like this is could be a convergence issue of the EM algorithm, but I don't know what one can do to overcome such problems. I've tried increasing the max_iter parameter to 20, but this doesn't seem to make much difference.

Would this be related to the quality of the data itself or to the fitting procedure?

It is not clear to me whether this is an issue with the data itself (or pre-processing), or with how I run scv.tl.recover_dynamics.

Starting with a scanpy-clustered dataset with the spliced/unspliced genes merged into an anndata object, my basic pipeline is just:

scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
scv.tl.recover_dynamics(adata, max_iter=20, n_jobs=4) # default: use_raw = False
scv.tl.velocity(adata, mode='dynamical', filter_genes = False) # default: use_raw = False, filter_genes = True
top_genes = adata.var['fit_likelihood'].sort_values(ascending=False).index
scv.pl.scatter(adata, color=cluster_name, basis=top_genes[:5], ncols=5, frameon=False)

Versions

```pytb scvelo==0.2.3 scanpy==1.8.1 anndata==0.7.6 loompy==3.0.6 numpy==1.20.1 scipy==1.6.0 matplotlib==3.4.3 sklearn==0.24.2 pandas==1.3.2 ```

[YitengDang]

I'm still looking for a solution to this problem. I read up a bit more on the EM algorithm, and noticed it can converge onto local maxima of the likelihood. I'm suspect that plots I'm including above seem might have this problem. Is there any way of preventing this, e.g. by changing or randomizing the initial conditions?