aristoteleo/dynamo-release

Inclusive model of expression dynamics with conventional or metabolic labeling based scRNA-seq / multiomics, vector field reconstruction and differential geometry analyses

dynamo-release's Issues

How to control the directionality?

Hi Xiaojie,

I tested dynamo on a time series data set. I get a complete reserve direction of the arrows since I have actual time points as a reference.

1. So i am wondering whether there is a way to control the directionality (such as reverse =T in monocle setting)?

2. Is there any quantitative way to measure dynamo is better than traditional RNA velocyto analysis? In here dynamo-notebooks/main_figs/Fig1C.ipynb , the visualization of the dynamo seems better ordered that the traditional RNA velocyto but sometimes it is kind of a little bit arbitrary.

Thanks in advance,
Peng

design a class hierarchy specific for dynamo

design a class hierarchy that includes molecule base class, then RNA, protein classes which inherit the base class where RNA can have labelling attribute or splicing attribute. We then have a single cell class that contains both RNA protein. refactor the entire code based on this strategy. Anndata will be just a data container for this class to input or ouput the content in this cluster

Adding the dataset for scNT-seq to sample_data

I was going through the dynamo tutorials and found it very helpful to play with the zebrafish dataset. Would it be possible to add the scNT-seq already preprocessed data to the sample_data in dynamo?

Thanks

An alternative way to show gene expression dynamics over time

Thanks for developing this powerful software to analyze the dynamics of cells from single-cell data.

We recently used dynamo to analyze the CART cell clone dynamics from single-cell time course data. It helped us a lot to elucidate accelerated and diverted cell differentiation dynamics. However, there are several issues that we could like to discuss to enable more people to use this powerful tool for more convenient analysis.

1. An alternative way to show gene expression dynamics over time
   Currently, we have gene expression dynamics over time like this:
   

Could we also offer an alternative way to show gene expression dynamics over time as a heatmap

Thus, we can compare the multiple genes in a uniform and formalized way.

Animate fate transition for PCA plot

Hi Xiaojie,

Thanks for this excellent tool. Is there a way to perform animate fate transition for PCA plot, not UMAP? Please see below, I tried, but it not work for me,

ax = dyn.pl.streamline_plot(adata, color='CellTypes', basis='pca', ax=ax, save_show_or_return='return')
instance = dyn.mv.StreamFuncAnim(adata=adata, basis='pca', color='CellTypes', ax=ax)
...........
raise ValueError('XA and XB must have the same number of columns '
ValueError: XA and XB must have the same number of columns (i.e. feature dimension.)

Thanks in advance.
Lin

Better dealing the relationship between layer/basis, velocity_key, vector field key, etc.

Currently, dynamo doesn't have nice convention for dealing with the relationship between layer/basis, velocity_key, vector field key, etc. By relationship, I mean each layer, each basis, should have corresponding velocity key, velocity field key, etc. A better strategy should be implemented.

add an utility function to identify bad genes with ntr goes up for a degradation experiment

If the ntr goes up for a gene in a degradation experiment, this may indicates some technical issues because the RNA should always decay for a chase experiment. Add a function to identify and filter those genes

make dynamo to support cell-wise alpha, beta and gamma estimation

transcription, splicing, and degradation rates should be gene, cell state-dependent. The current estimation framework for labeling data can be conveniently extended to estimate cell-wise alpha, beta and gamma rates (although most of beta and gamma are more or less constant). It is also possible to estimate cell-wise alpha, beta and gamma rates for conventional scRNA-seq data.

dynamo-notebooks/NASC_seq/NASC_seq_analysis.ipynb SparseVFC error

Hello, dynamo team,

When I am using the latest dynamo to run your jupyter notebook, it output the following errors for the dyn.tl.VectorField function. The notebook seems to be not compatible with the latest dynamo.

UnboundLocalError: local variable 'P' referenced before assignment.

Could we output the dynamics as a network with growth and transition rates?

Thanks for developing this powerful software to analyze the dynamics of cells from single-cell data.

We recently used dynamo to analyze the CART cell clone dynamics from single-cell time course data. It helped us a lot to elucidate accelerated and diverted cell differentiation dynamics. However, there are several issues that we could like to discuss to enable more people to use this powerful tool for more convenient analysis.

Currently, the RNA velocity output is mostly a visualization. When we want to quantify the dynamics with some numbers, could we output the dynamics as a network with growth and transition rates? Such as:

Update the the contributing.md file

Dynamo welcomes contribution (documentation improvement, bug fix, new features / test, organizational improvement, code readability, etc) from the community. The contribution guidance needs to be updated to make this process more welcoming and easy to follow.

cleanup code to ensure all code is compatible with black and flake

Currently, there are many places in our code that don't follow the black and flake convention. We need to update our code to ensure all codes follow those conventions.

add a new function to plot the nb model fitting

add a new function to plot the NB model fitting

Can't use dyn.pl.topography with PCA vector fields dim >2

Hi,

Thanks for this great tool, I've really enjoyed exploring the new functionalities (especially fate animation!). I've noticed an issue regarding trying to plot topographies of vector fields in PCA spaces of more than 2 dimensions.

If I run through the pipeline delineated in the Zebrafish tutorial with UMAP as my choice of basis, it works fine whether i specify 2 or >2 dimensions in dyn.pp.recipe_monocle() and dyn.tl.reduceDimension(). If I exclusively use PCA as my basis, and specify exactly 2 dimensions, it also works fine. However, if I try to use a PCA basis with >2 dimensions, dyn.vf.VectorField() works fine but when I try to run dyn.pl.topography() I get the following error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-148-3142c99fe9de> in <module>
----> 1 dyn.pl.topography(dyn_data, basis='pca', color='cell_subtype_unique', terms=['streamline'],  x=0, y=1)

~/.local/lib/python3.8/site-packages/dynamo/plot/topography.py in topography(adata, basis, x, y, color, layer, highlights, labels, values, theme, cmap, color_key, color_key_cmap, background, ncols, pointsize, figsize, show_legend, use_smoothed, xlim, ylim, t, terms, init_cells, init_states, quiver_source, fate, approx, quiver_size, quiver_length, density, linewidth, streamline_color, streamline_alpha, color_start_points, markersize, marker_cmap, save_show_or_return, save_kwargs, aggregate, show_arrowed_spines, ax, sort, frontier, s_kwargs_dict, q_kwargs_dict, **streamline_kwargs_dict)

~/.local/lib/python3.8/site-packages/dynamo/vectorfield/topography.py in topography(adata, basis, layer, X, dims, n, VecFld)

~/.local/lib/python3.8/site-packages/dynamo/vectorfield/topography.py in find_fixed_points_by_sampling(self, n, x_range, y_range, lhs, tol_redundant)

~/.local/lib/python3.8/site-packages/dynamo/vectorfield/topography.py in find_fixed_points(X0, func_vf, tol_redundant, full_output)

/camp/lab/briscoej/working/Rory/.conda/my_envs/py_lab/lib/python3.8/site-packages/scipy/optimize/minpack.py in fsolve(func, x0, args, fprime, full_output, col_deriv, xtol, maxfev, band, epsfcn, factor, diag)
    145                'diag': diag}
    146 
--> 147     res = _root_hybr(func, x0, args, jac=fprime, **options)
    148     if full_output:
    149         x = res['x']

/camp/lab/briscoej/working/Rory/.conda/my_envs/py_lab/lib/python3.8/site-packages/scipy/optimize/minpack.py in _root_hybr(func, x0, args, jac, col_deriv, xtol, maxfev, band, eps, factor, diag, **unknown_options)
    211     if not isinstance(args, tuple):
    212         args = (args,)
--> 213     shape, dtype = _check_func('fsolve', 'func', func, x0, args, n, (n,))
    214     if epsfcn is None:
    215         epsfcn = finfo(dtype).eps

/camp/lab/briscoej/working/Rory/.conda/my_envs/py_lab/lib/python3.8/site-packages/scipy/optimize/minpack.py in _check_func(checker, argname, thefunc, x0, args, numinputs, output_shape)
     24 def _check_func(checker, argname, thefunc, x0, args, numinputs,
     25                 output_shape=None):
---> 26     res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
     27     if (output_shape is not None) and (shape(res) != output_shape):
     28         if (output_shape[0] != 1):

~/.local/lib/python3.8/site-packages/dynamo/vectorfield/utils_vecCalc.py in <lambda>(x)

~/.local/lib/python3.8/site-packages/dynamo/tools/utils.py in timed(*args, **kw)

~/.local/lib/python3.8/site-packages/dynamo/vectorfield/utils_vecCalc.py in vector_field_function(x, vf_dict, dim, kernel, **kernel_kwargs)

~/.local/lib/python3.8/site-packages/dynamo/tools/utils.py in timed(*args, **kw)

~/.local/lib/python3.8/site-packages/dynamo/vectorfield/utils_vecCalc.py in con_K(x, y, beta, method, return_d)

/camp/lab/briscoej/working/Rory/.conda/my_envs/py_lab/lib/python3.8/site-packages/scipy/spatial/distance.py in cdist(XA, XB, metric, *args, **kwargs)
   2719         raise ValueError('XB must be a 2-dimensional array.')
   2720     if s[1] != sB[1]:
-> 2721         raise ValueError('XA and XB must have the same number of columns '
   2722                          '(i.e. feature dimension.)')
   2723 

ValueError: XA and XB must have the same number of columns (i.e. feature dimension.)

It's hard for me to parse exactly where things are breaking down - if I might be doing something wrong let me know!

Cheers,
Rory

P.s.

If you need more details about the code I ran let me know, here's my dependency info:

package	matplotlib	numba	numpy	pandas	python-igraph	scikit-learn	scipy	seaborn	setuptools	statsmodels	tqdm	hdbscan	cvxopt	trimap	numdifftools	colorcet	anndata	dynamo-release	loompy	pynndescent	umap-learn
version	3.2.0	0.48.0	1.18.1	1.0.2	0.8.0	0.22.2.post1	1.4.1	0.10.0	46.0.0.post20200311	0.11.1	4.43.0	0.8.26	1.2.4	1.4.3.dev1	0.9.39	2.0.2	0.7.3	0.95.1	3.0.6	0.4.8	0.4.3

Failed to install package, because of statsmodels

Hi,
Thanks for providing this exciting single-cell analysis tool. I have some problems while installing dynamo. Now I am using Python 3.7.3 with miniconda. I am sure I have installed the package statsmodels (version 0.10.1), which can be loaded successfully in the python session. But I always got the error message while installing the dynamo " import statsmodels.api as sm
ModuleNotFoundError: No module named 'statsmodels'". Anything could you suggest? For example, is there any specific version of statsmodels needed?
Thank you!!

Matplotlib issue

Hi,

We were testing the package on a cluster with python 3.8 and when it comes to

dyn.pl.streamline_plot(adata, save_show_or_return='save',save_kwargs={"prefix":'dynamo.streamline',"dpi":1000, "ext":'png'},contour=True)

it pops up "In:" without a way to exit.

I guess this is an issue with matplotlib because when we try to set

dyn.configuration.set_figure_params('dynamo', background='white')

The same issue pops out again.

The structure of adata is:

adata.obs_keys
<bound method AnnData.obs_keys of AnnData object with n_obs × n_vars = 7525 × 60662
obs: 'TotalUMIs', 'orig_ident', 'batch', 'nGenes', 'nCounts', 'pMito', 'use_for_pca', 'spliced_Size_Factor', 'initial_spliced_cell_size', 'unspliced_Size_Factor', 'initial_unspliced_cell_size', 'Size_Factor', 'initial_cell_size', 'ntr', 'cell_cycle_phase', 'jaccard_velocity_confidence', 'umap_ddhodge_sampled', 'umap_ddhodge_div', 'umap_potential'
var: 'Accession', 'AccessionVersion', 'Aliases', 'CcdsID', 'Chromosome', 'ChromosomeEnd', 'ChromosomeStart', 'CosmicID', 'DnaBindingDomain', 'FullName', 'GeneType', 'HgncID', 'IsTF', 'Location', 'LocationSortable', 'LocusGroup', 'LocusType', 'MgdID', 'MirBaseID', 'OmimID', 'PubmedID', 'RefseqID', 'RgdID', 'UcscID', 'UniprotID', 'VegaID', 'pass_basic_filter', 'score', 'log_cv', 'log_m', 'use_for_pca', 'ntr', 'beta', 'gamma', 'half_life', 'alpha_b', 'alpha_r2', 'gamma_b', 'gamma_r2', 'gamma_logLL', 'delta_b', 'delta_r2', 'uu0', 'ul0', 'su0', 'sl0', 'U0', 'S0', 'total0', 'use_for_dynamics', 'use_for_transition'
uns: 'velocyto_SVR', 'pp_norm_method', 'PCs', 'explained_variance_ratio_', 'pca_fit', 'feature_selection', 'dynamics', 'neighbors', 'umap_fit', 'grid_velocity_umap', 'VecFld_umap', 'VecFld'
obsm: 'X_pca', 'X', 'cell_cycle_scores', 'X_umap', 'velocity_umap', 'velocity_umap_SparseVFC', 'X_umap_SparseVFC'
layers: 'matrix', 'spliced', 'unspliced', 'X_spliced', 'X_unspliced', 'M_u', 'M_uu', 'M_s', 'M_us', 'M_ss', 'velocity_S'
obsp: 'moments_con', 'connectivities', 'distances', 'pearson_transition_matrix'>

We were using it like:

import dynamo as dyn
import anndata
import pandas as pd
import numpy as np
import matplotlib
import matplotlib.pyplot as plt

Missing arguments for dyn.tl.dynamics

Python 3.8.5 (default, Sep  4 2020, 07:30:14) 
[GCC 7.3.0] :: Anaconda, Inc. on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import dynamo as dyn
>>> adata = dyn.read('/home/tjcooper/Projects/MDSC_Project/mouseBM.h5ad')
>>> dyn.pp.recipe_monocle(adata)
/home/tjcooper/dynamo-release/dynamo/preprocessing/utils.py:439: UserWarning: `total_szfactor` is not `None` and it is not in adata object.
  warnings.warn("`total_szfactor` is not `None` and it is not in adata object.")
AnnData object with n_obs × n_vars = 2865 × 14836
    obs: 'orig.ident', 'nCount_spliced', 'nFeature_spliced', 'nCount_unspliced', 'nFeature_unspliced', 'nCount_ambiguous', 'nFeature_ambiguous', 'nCount_RNA', 'nFeature_RNA', 'nCount_SCT', 'nFeature_SCT', 'seurat_clusters', 'nGenes', 'nCounts', 'pMito', 'use_for_pca', 'Size_Factor', 'initial_cell_size', 'unspliced_Size_Factor', 'initial_unspliced_cell_size', 'spliced_Size_Factor', 'initial_spliced_cell_size', 'ntr', 'cell_cycle_phase'
    var: 'features', 'ambiguous_features', 'spliced_features', 'unspliced_features', 'nCells', 'nCounts', 'pass_basic_filter', 'log_cv', 'score', 'log_m', 'use_for_pca', 'ntr'
    uns: 'pp', 'velocyto_SVR', 'PCs', 'explained_variance_ratio_', 'pca_fit', 'feature_selection'
    obsm: 'X_umap', 'X_pca', 'X', 'cell_cycle_scores'
    layers: 'ambiguous', 'spliced', 'unspliced', 'X_unspliced', 'X_spliced'
>>> dyn.tl.dynamics(adata)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/tjcooper/dynamo-release/dynamo/tools/dynamics.py", line 303, in dynamics
    moments(adata, genes=valid_bools, group=tkey)
  File "/home/tjcooper/dynamo-release/dynamo/tools/moments.py", line 99, in moments
    kNN, knn_indices, knn_dists, _ = umap_conn_indices_dist_embedding(
  File "/home/tjcooper/dynamo-release/dynamo/tools/connectivity.py", line 227, in umap_conn_indices_dist_embedding
    embedding_ = simplicial_set_embedding(
TypeError: simplicial_set_embedding() missing 3 required positional arguments: 'densmap', 'densmap_kwds', and 'output_dens'

I'm using the PyPi release of Dynamo, and receiving the above error message when attempting to run it on this file:

https://e.pcloud.link/publink/show?code=XZuibkZCBtw4neJ5Q8wpFqcCoHUKjBx2dTX

Preprocessing was originally performed with Seurat and the h5ad file was created via SeuratDisk.

I've also tried reading in the original, unprocessed .loom file via scanpy.read_loom but encounter an identical error:

https://e.pcloud.link/publink/show?code=XZvvbkZfK2Q2yvHnhS8YA0NCtXDWLzqQwh7

What could be causing this?

About data in scNT-seq

Dear Dr Qiu,

Thank you very much for this great package. I found you use tow data sets in scNT-seq. One is splicing dataset. Another is labeling dataset. Splicing dataset can get the spliced and unspliced count. Labeling dataset can get new counts and total counts. I knew the labeling dataset is the result of scNT-seq. But I really want to know how can I get the splicing dataset. Is it another experiment?

An alternative way to show gene expression dynamics over time.

Thanks for developing this powerful software to analyze the dynamics of cells from single-cell data.

We recently used dynamo to analyze the CART cell clone dynamics from single-cell time course data. It helped us a lot to elucidate accelerated and diverted cell differentiation dynamics. However, there are several issues that we could like to discuss to enable more people to use this powerful tool for more convenient analysis.

1. An alternative way to show gene expression dynamics over time
   Currently, we have gene expression dynamics over time like this:
   

Could we also offer an alternative way to show gene expression dynamics over time as a heatmap

Thus, we can compare the multiple genes in a uniform and formalized way.

Vectorfield needs more online documentation

Thanks for making API page online.
Currently there is no online documentation about vectorfield object after running vf.VectorField, though source code comments explain this result object stored in adata well. Meanwhile, I would like to apply an existing vectorfield to some new data point to do some theoretical stuff. For now my pipeline calls vectorfield.func(X). It is handy if there is a utility function apply_vectorfield(adata_new, vectorfield) or apply_vectorfield(X, vectorfield).

Meanwhile, is topography of vectorfield only supporting 2D? If I use PCA basis, then there is some error raised.

Use logging in dyn.vf.scVectorField.py instead of verbose level

We would like to replace the printing of running information in vector field reconstruction via verbose with the logging system from dynamo:
def SparseVFC( X: np.ndarray, Y: np.ndarray, Grid: np.ndarray, M: int = 100, a: float = 5, beta: float = None, ecr: float = 1e-5, gamma: float = 0.9, lambda_: float = 3, minP: float = 1e-5, MaxIter: int = 500, theta: float = 0.75, div_cur_free_kernels: bool = False, velocity_based_sampling: bool = True, sigma: float = 0.8, eta: float = 0.5, seed=0, lstsq_method: str = "drouin", verbose: int = 1,)

Error with simplicial set embedding

Hi,

Thanks for this great package - I understand that it is still under development so apologies if this error is something that is already known about.

I'm running dynamo version '0.0+11bf3a9' and Python 3.7.3. For a very simple test case based on sample data, I run:

import dynamo as dyn
adata = dyn.sample_data.BM()
dyn.pp.recipe_monocle(adata, n_top_genes=1000)

dyn.tl.dynamics(adata, filter_gene_mode="no", use_smoothed=True,
    experiment_type='conventional', assumption_mRNA='ss', model='deterministic',
    est_method='linear_regression', log_unnormalized=True)

The error I get is 'AttributeError: "tuple' object has no attribute 'tocoo' ". This occurs in line 201 of connectivity.py, in umap_conn_indices_dist_embedding(), when umap's simplicial_set_embedding() is called.

It appears that umap's simplicial_set_embedding() doesn't just output a scipy sparse matrix, but outputs a tuple, of which the first object is the desired matrix object (my umap version is '0.4.3').

I thought it would be worth raising this issue as it could be a simple solution, or something simple going wrong for me!

Cheers!

Cannot install package

Hi dynamo group,

I used bellow command to install however it didn't work for me. Is something wrong with my set-up?
pip install git+https://github.com:aristoteleo/dynamo-release

Collecting git+https://github.com:aristoteleo/dynamo-release
Cloning https://github.com:aristoteleo/dynamo-release to /tmp/pip-req-build-i1smn3l5
Running command git clone -q https://github.com:aristoteleo/dynamo-release /tmp/pip-req-build-i1smn3l5
fatal: unable to access 'https://github.com:aristoteleo/dynamo-release/': Port number ended with 'a'
ERROR: Command "git clone -q https://github.com:aristoteleo/dynamo-release /tmp/pip-req-build-i1smn3l5" failed with error code 128 in None

Best,
Peng

Is it possible to choose a subset of cells to do a streamline plot?

Thanks for developing this powerful software to analyze the dynamics of cells from single-cell data.

We recently used dynamo to analyze the CART cell clone dynamics from single-cell time course data. It helped us a lot to elucidate accelerated and diverted cell differentiation dynamics. However, there are several issues that we could like to discuss to enable more people to use this powerful tool for more convenient analysis.

1. Currently streamline plot the whole cell population, is it possible to choose a subset of cells to do a streamline plot? For example, we have several clones of cells from our lineage tracing, could we plot a streamline plot from each clone?
   

2. Could we plot a transition-growth network for each clone as well?
   

[related to paper]What is cscRNA-seq

I couldn't find the definition of cscRNA-seq. Do you know what it is?

Add FAQ page in readthedocs

There are a lot of frequently asked questions (FAQs) regarding dynamo from the community. We should add a FAQ page in our docs for the benefits of our great community!

From Seurat to Dynamo

Hello,

I would be interested in trying Dynamo, however, all the data pre-processing of my data has been done in seurat and scater. Would it be feasible to start with Dynamo from a pre-processed dataset? (I also have a loom file as output of velocyto)

Many thanks

module 'dynamo.tl' has no attribute 'estimation'

Hi @Xiaojieqiu ,
I followed the dynamo-notebooks/NASC_seq/NASC_seq_analysis.ipynb ( In [3] ).

The command "dyn.tl.estimation" returned "module 'dynamo.tl' has no attribute 'estimation'"
Also "dyn.tl.velocity" returned "module 'dynamo.tl' has no attribute 'velocity'"

Did I install dynamo unsucseccfully or missed some commands?

Thanks,
Raulee

Error running dyn.tl.moments()

Good day!

Thanks for developing this great tool! I am trying to analyze some of my data but had a problem when running dyn.tl.moments(). I have exported my Seurat object as loom file and using scvelo, I have read and merged the new anndata object with the splicing info.

AnnData object with n_obs × n_vars = 7116 × 25639
    obs: 'AC3', 'ClusterID', 'ClusterName', 'OC2', 'OPC_shared4', 'OPC_variable5', 'RNA_snn_res_0_3', 'RNA_snn_res_0_7', 'cellcycle1', 'nCount_RNA', 'nFeature_RNA', 'orig_ident', 'seurat_clusters', 'sample_batch', 'Clusters', '_X', '_Y', 'initial_size_spliced', 'initial_size_unspliced', 'initial_size'
    var: 'Selected', 'vst_mean', 'vst_variable', 'vst_variance', 'vst_variance_expected', 'vst_variance_standardized', 'Accession', 'Chromosome', 'End', 'Start', 'Strand'
    obsm: 'mnn_cell_embeddings', 'pca_cell_embeddings', 'umap_cell_embeddings'
    varm: 'pca_feature_loadings'
    layers: 'norm_data', 'scale_data', 'ambiguous', 'matrix', 'spliced', 'unspliced'

Then, I applied the dyn.pp.recipe_monocle(adata):

AnnData object with n_obs × n_vars = 7116 × 25639
    obs: 'AC3', 'ClusterID', 'ClusterName', 'OC2', 'OPC_shared4', 'OPC_variable5', 'RNA_snn_res_0_7', 'cellcycle1', 'nCount_RNA', 'nFeature_RNA', 'orig_ident', 'seurat_clusters', 'sample_batch', 'Clusters', '_X', '_Y', 'initial_size_spliced', 'initial_size_unspliced', 'initial_size', 'nGenes', 'nCounts', 'pMito', 'use_for_pca', 'scale_data_Size_Factor', 'initial_scale_data_cell_size', 'spliced_Size_Factor', 'initial_spliced_cell_size', 'norm_data_Size_Factor', 'initial_norm_data_cell_size', 'unspliced_Size_Factor', 'initial_unspliced_cell_size', 'Size_Factor', 'initial_cell_size', 'ntr', 'cell_cycle_phase'
    var: 'Selected', 'vst_mean', 'vst_variable', 'vst_variance', 'vst_variance_expected', 'vst_variance_standardized', 'Accession', 'Chromosome', 'End', 'Start', 'Strand', 'pass_basic_filter', 'log_cv', 'score', 'log_m', 'use_for_pca', 'ntr'
    uns: 'velocyto_SVR', 'pp_norm_method', 'PCs', 'explained_variance_ratio_', 'pca_fit', 'feature_selection'
    obsm: 'mnn_cell_embeddings', 'pca_cell_embeddings', 'umap_cell_embeddings', 'X_pca', 'X', 'cell_cycle_scores'
    varm: 'pca_feature_loadings'
    layers: 'norm_data', 'scale_data', 'ambiguous', 'matrix', 'spliced', 'unspliced', 'X_scale_data', 'X_spliced', 'X_norm_data', 'X_unspliced'

But when trying to run either first dyn.tl.moments(adata) or dyn.tl.dynamics(adata) directly, I get the error below:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-9-47e88b9d5c17> in <module>
----> 1 dyn.tl.moments(adata)

/hpc/pmc_stunnenberg/cruiz/miniconda3/envs/python3/lib/python3.8/site-packages/dynamo/tools/moments.py in moments(adata, genes, group, use_gaussian_kernel, normalize, use_mnn, layers, n_pca_components, n_neighbors)
    148             layer_y = adata.layers[layer2].copy()
    149 
--> 150             layer_y_group = np.where([layer2 in x for x in
    151                                       [only_splicing, only_labeling, splicing_and_labeling]])[0][0]
    152             # don't calculate 2 moments among uu, ul, su, sl -

IndexError: index 0 is out of bounds for axis 0 with size 0

Do you know what the problem might be? Thanks in advance for your help!

ensure all transition matrix are row normalized

the transition matrix from the correlation kernel is row normalized while transition matrices from markov chain classes assume column normalization. This should be unified.

concatenate and layers for scEU-seq data

When we tried to import the scEU-seq data from the two libraries, labelled and unlabelled, using anndata.concat(),
we have encountered some problems reported that InvalidIndexError: Reindexing only valid with uniquely valued Index objects.
So, we seek to direct output the merged libraries via using one library as basis and borrowing layers from another, and it goes like this

merged=unlablled.copy()
merged.layers['ls']=labelled.layers['spliced']
merged.layers['lu']=labelled.layers['unspliced']
merged.layers['us']=unlabelled.layers['spliced']
merged.layers['uu']=unlabelled.layers['unspliced']
del merged.layers['spliced']
del merged.layers['unspliced']
del merged.layers['matrix']
del merged.layers['ambiguous']
merged.layers['new']=merged.layers['ls']+merged.layers['lu']
merged.layers['total']=merged.layers['ls']+merged.layers['lu']+merged.layers['us']+merged.layers['uu']
del merged.layers['ls']
del merged.layers['lu']
del merged.layers['us']
del merged.layers['uu']
merged

but when we examined the merged files with dyn.pl.basic_stats(), it reports only the unlabelled data
So, we are wondering an appropriate way to merge the two libraries and import them as dynamo input

Import error: hdbscan initialization error on Python >= 3.7.9 and Numpy <= 1.18.1

'ValueError: numpy.ndarray size changed, may indicate binary incompatibility. Expected 88 from C header, got 80 from PyObject'

Strongly recommend updating the min Numpy version requirement to 1.20.0.

ZeroDivisionError: division by zero

We need to ensure all functions in our package will have special consideration when the denominator is 0, especially for the functions in kinetic parameter estimation.

Possible bug in find_fixed_points() (+ tqdm progress bars)

Hey,

I've encounter a strange error a couple times when running dyn.vf.VectorField(). It doesn't happen every time, and is subject to the particular data and dimensionality of the vector field being learnt. As such, it is pretty hard for me to give thoroughly reproducible code, but it seems to me that the issue happens only in certain contexts depending on the output of remove_redundant_points() within find_fixed_points().

Sometime when I run something like:
vf = dyn.vf.VectorField(adata, basis='pca', return_vf_object=True, pot_curl_div=True, dims=[0,1])

I get the following error message:

s)
    207             data_2d = self.adata.copy()
    208             print('Low-dimension:')
--> 209             vf_2d = dyn.vf.VectorField(data_2d, basis=self.basis, return_vf_object=True, pot_curl_div=True, dims=[0,1])
    210             self.vf_dict['{0}_2D'.format(name)] = (data_2d, vf_2d)
    211             if plot:

~/.local/lib/python3.8/site-packages/dynamo_release-0.95.2-py3.8.egg/dynamo/vectorfield/topography.py in VectorField(adata, basis, layer, dims, genes, normalize, grid_velocity, grid_num, velocity_key, method, return_vf_object, map_topography, pot_curl_div, cores, **kwargs)
    701             warnings.simplefilter("ignore")
    702 
--> 703             adata = topography(
    704                 adata, basis=basis, X=X, layer=layer, dims=[0, 1], VecFld=vf_dict['VecFld'], **tp_kwargs
    705             )

~/.local/lib/python3.8/site-packages/dynamo_release-0.95.2-py3.8.egg/dynamo/vectorfield/topography.py in topography(adata, basis, layer, X, dims, n, VecFld)
    499     vecfld.find_fixed_points_by_sampling(n, xlim, ylim)
    500     if vecfld.get_num_fixed_points() > 0:
--> 501         vecfld.compute_nullclines(xlim, ylim, find_new_fixed_points=True)
    502     # sep = compute_separatrices(vecfld.Xss.get_X(), vecfld.Xss.get_J(), vecfld.func, xlim, ylim)
    503     #

~/.local/lib/python3.8/site-packages/dynamo_release-0.95.2-py3.8.egg/dynamo/vectorfield/topography.py in compute_nullclines(self, x_range, y_range, find_new_fixed_points, tol_redundant)
    439         if find_new_fixed_points:
    440             sample_interval = ds * 10
--> 441             X, J = find_fixed_points_nullcline(
    442                 self.func, self.NCx, self.NCy, sample_interval, tol_redundant, True
    443             )

~/.local/lib/python3.8/site-packages/dynamo_release-0.95.2-py3.8.egg/dynamo/vectorfield/topography.py in find_fixed_points_nullcline(func, NCx, NCy, sample_interval, tol_redundant, full_output)
    227     int_P = np.array(int_P)
    228     if full_output:
--> 229         P, J, _ = find_fixed_points(int_P, func, tol_redundant, full_output=True)
    230         return P, J
    231     else:

~/.local/lib/python3.8/site-packages/dynamo_release-0.95.2-py3.8.egg/dynamo/vectorfield/topography.py in find_fixed_points(X0, func_vf, tol_redundant, full_output)
     61         if full_output:
     62             X, discard = remove_redundant_points(X, tol_redundant, output_discard=True)
---> 63             J = J[~discard]
     64             fval = fval[~discard]
     65         else:

IndexError: boolean index did not match indexed array along dimension 0; dimension is 0 but corresponding boolean dimension is 1

It seems like an error with numpy dimensions but I'm not sure exactly what's causing it. Sorry I couldn't be more helpful this time. Let me know if you know what's causing this.

Also, since I'm posting I thought I'd drop an additional little request: I see that with fate prediction and vector field learning, tqdm progress bars are used. This verbosity is actually a bit undesirable for my usage case, so it would be excellent to have a user argument to disable progress bars/verbosity.

Thanks again for this fantastic package, I've really enjoyed using it so far.

Cheers,
Rory

How to visualize development trajectory with deeper channel?

Thanks for developing this powerful software to analyze the dynamics of cells from single-cell data.

We recently used dynamo to analyze the CART cell clone dynamics from single-cell time course data. It helped us a lot to elucidate accelerated and diverted cell differentiation dynamics. However, there are several issues that we could like to discuss to enable more people to use this powerful tool for more convenient analysis.

1. How to visualize development trajectory with a deeper channel?

From our experiment data, we know that one sample of cells differentiate 3 days earlier and cells are more committed. It seems that the differentiation channel is deeper and steeper. Could we have a better way to visualize and distinguish these two cases?

2. Usually top 30 genes will determine the RNA velocity. Could we identify the genes which control the lineage decision and contribute to certain cell lineage most? It is a little bit like the toggle switch genes near the branching point which are identified by Monocle.

Inverse transformation of fate predictions giving errors

Hi Xiaojie,

Exploring dynamo's fate prediction functions, I find that I am unable to run dyn.pd.fate() with inverse_transform=True.

Working with a PCA space of 30 dims, I run the following command:

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

and I get the following error message:

integration with ivp solver: 100%|██████████| 271/271 [00:13<00:00, 20.06it/s]
uniformly sampling points along a trajectory: 100%|██████████| 271/271 [00:01<00:00, 162.27it/s]
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-15-35417f0465e8> in <module>
----> 1 dyn.pd.fate(adata, basis='pca', init_cells=progenitor, interpolation_num=100,  direction='forward',
      2    inverse_transform=True, average=False, cores=1)

/camp/lab/briscoej/working/Rory/.conda/my_envs/py_lab/lib/python3.8/site-packages/dynamo/prediction/fate.py in fate(adata, init_cells, init_states, basis, layer, dims, genes, t_end, direction, interpolation_num, average, sampling, VecFld_true, inverse_transform, scale, cores, **kwargs)
    145     high_prediction = None
    146     if basis == "pca" and inverse_transform:
--> 147         high_prediction = adata.uns["pca_fit"].inverse_transform(prediction)
    148         if adata.var.use_for_dynamics.sum() == high_prediction.shape[1]:
    149             valid_genes = adata.var_names[adata.var.use_for_dynamics]

/camp/lab/briscoej/working/Rory/.conda/my_envs/py_lab/lib/python3.8/site-packages/sklearn/decomposition/_base.py in inverse_transform(self, X)
    157                             self.components_) + self.mean_
    158         else:
--> 159             return np.dot(X, self.components_) + self.mean_

<__array_function__ internals> in dot(*args, **kwargs)

ValueError: could not broadcast input array from shape (30,99) into shape (30)

It looks like prediction is of shape (n_components, n_datapoints) but for inverse_transform() it should be of shape (n_datapoints, n_components). Hopefully this should be a quick fix!

Thanks again for a great package.

Best,
Rory

Here are my dependencies:

package 	anndata	dynamo-release	matplotlib	numba	numpy	pandas	python-igraph	scikit-learn	scipy	seaborn	setuptools	statsmodels	tqdm	hdbscan	cvxopt	trimap	numdifftools	colorcet	loompy	pynndescent	umap-learn
version	0.7.4	0.95.2	3.2.0	0.48.0	1.18.1	1.0.2	0.8.0	0.22.2.post1	1.4.1	0.10.0	46.0.0.post20200311	0.11.1	4.43.0	0.8.26	1.2.4	1.4.3.dev1	0.9.39	2.0.2	3.0.6	0.4.8	0.4.3

P.s. in case helpful: after loading my own data set the code I run is:

dyn.pp.recipe_monocle(adata)
dyn.tl.dynamics(adata)
dyn.tl.moments(adata)
dyn.tl.reduceDimension(adata, basis='pca')
dyn.tl.cell_velocities(adata, basis='pca')
dyn.vf.VectorField(adata, basis='pca')
progenitor = adata[adata.obs['cell_type']=='Progenitor'].obs_names
dyn.pd.fate(adata, basis='pca', init_cells=progenitor, interpolation_num=100,  direction='forward',  inverse_transform=True, average=False, cores=1)

add an util function to detect whether the nearest neighbor graph is broken

After subsetting adata object, the nearest neighbor graph need to be recomputed. Otherwise ddhoge (and other functions potentially) may break. We need a function to detect the completeness of the nearest neighbor graph and then recompute them for all functions that uses the nearest neighbor graph.

add a new file (clusters.py in tools) to use cdlib for various community clustering

We would like to improve our support for cell clustering and annotation. To beginning with, let us start with using cdlib package to support various clustering algorithms in dynamo.

The interface should something like: dyn.tl.clusters(adata, method=method)

method can be: infomap, louvain, leiden, hdbscan, kmeans,

readthedocs: https://cdlib.readthedocs.io/en/latest/
notebook: https://colab.research.google.com/github/GiulioRossetti/cdlib/blob/master/docs/CDlib.ipynb#scrollTo=Oj_-Xhh_4bps

Is there a docker version for dynamo?

Hello,
I find dynamo a great tool for single-cell analysis tool.
However, I found it hard to install this package because the network of my cluster is bad.
Is there possible provide a docker image of dynamo?

Thanks!!!

How to deal with the error in phase plot?

When we are plotting a phage portrait, there is an error pop up. How to solve it?

Harmonize anndata layer nomenclature in a pulse - chase experiment

Hi Xiaojie,

Currently the input to dyn.tl.recipe_deg_data requires the layer with labeled counts to be named 'new'. To better reflect the nature of those counts in a pulse chase experiment, it could be good to use 'labeled' as the name for that layer instead.

Thanks,
Jorge

create a public slack channel

we can follow the instruction from this page to create a public slack channel that everyone can sign up.

I tried to add an invitation link to the readme.md before but this link will expire in one month.

Error when running dyn.tl.dynamics

I get the following error when running dyn.tl.dynamics on the DentateGyrus test data set:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/vfranke/bin/Software/miniconda3/envs/scvelo/lib/python3.8/site-packages/dynamo/tools/dynamics.py", line 282, in dynamics
    moments(adata, genes=valid_bools, group=tkey)
  File "/home/vfranke/bin/Software/miniconda3/envs/scvelo/lib/python3.8/site-packages/dynamo/tools/moments.py", line 95, in moments
    kNN, knn_indices, knn_dists, _ = umap_conn_indices_dist_embedding(
  File "/home/vfranke/bin/Software/miniconda3/envs/scvelo/lib/python3.8/site-packages/dynamo/tools/connectivity.py", line 227, in umap_conn_indices_dist_embedding
    embedding_ = simplicial_set_embedding(
TypeError: simplicial_set_embedding() missing 3 required positional arguments: 'densmap', 'densmap_kwds', and 'output_dens'

package versions:

scvelo 0.2.3
dynamo-release 0.95.2
Python 3.8.8

CLT for quantifying mutation -- need update:

1. filter reads that have mutations and send them to velocyto or the recent bus tool.
2. create the u_u, u_l, s_u, s_l matrix
3. process CITE-seq, REAP-seq data also

update our contributing file

We need to update our contributing guidelines so that they can guide other developers.

error when install this program

when I use the code "python3 -m pip install git+https://github.com:aristoteleo/dynamo-release" to install this software, but the error was rraised as follow.
Collecting git+https://github.com:aristoteleo/dynamo-release
Cloning https://github.com:aristoteleo/dynamo-release to /tmp/pip-req-build-3yrf835a
Running command git clone -q https://github.com:aristoteleo/dynamo-release /tmp/pip-req-build-3yrf835a
fatal: unable to access 'https://github.com:aristoteleo/dynamo-release/': URL using bad/illegal format or missing URL
ERROR: Command errored out with exit status 128: git clone -q https://github.com:aristoteleo/dynamo-release /tmp/pip-req-build-3yrf835a Check the logs for full command output.
So, what your suggestion for this error?

Some problems when I save a PDF plot

I get some problems when I tried to save a plot to a PDF file. The plot looks good in the notebook, but I cannot save it to a PDF.

Saving figure to ./plot/celltype.pdf...

ValueError Traceback (most recent call last)
in
1 dyn.pl.streamline_plot(adata, color=['celltype'], basis='umap', show_legend='on data', show_arrowed_spines=True, pointsize = 0.3, theme = "darkred", save_show_or_return = "return")
----> 2 dyn.pl.save_fig("./plot/celltype")
3
4
~/.local/lib/python3.7/site-packages/dynamo_release-0.95.2-py3.7.egg/dynamo/plot/utils.py in save_fig(path, prefix, dpi, ext, transparent, close, verbose)
1251
1252 # Actually save the figure
-> 1253 plt.savefig(savepath, dpi=300 if dpi is None else dpi, transparent=transparent, format=ext, bbox_inches='tight')
1254
1255 # Close it
~/miniconda3/lib/python3.7/site-packages/matplotlib/pyplot.py in savefig(*args, **kwargs)
857 def savefig(*args, **kwargs):
858 fig = gcf()
--> 859 res = fig.savefig(*args, **kwargs)
860 fig.canvas.draw_idle() # need this if 'transparent=True' to reset colors
861 return res
~/miniconda3/lib/python3.7/site-packages/matplotlib/figure.py in savefig(self, fname, transparent, **kwargs)
2309 patch.set_edgecolor('none')
2310
-> 2311 self.canvas.print_figure(fname, **kwargs)
2312
2313 if transparent:
~/miniconda3/lib/python3.7/site-packages/matplotlib/backend_bases.py in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
2215 orientation=orientation,
2216 bbox_inches_restore=_bbox_inches_restore,
-> 2217 **kwargs)
2218 finally:
2219 if bbox_inches and restore_bbox:
~/miniconda3/lib/python3.7/site-packages/matplotlib/backend_bases.py in wrapper(*args, **kwargs)
1637 kwargs.pop(arg)
1638
-> 1639 return func(*args, **kwargs)
1640
1641 return wrapper
~/miniconda3/lib/python3.7/site-packages/matplotlib/backends/backend_pdf.py in print_pdf(self, filename, dpi, bbox_inches_restore, metadata)
2591 RendererPdf(file, dpi, height, width),
2592 bbox_inches_restore=bbox_inches_restore)
-> 2593 self.figure.draw(renderer)
2594 renderer.finalize()
2595 if not isinstance(filename, PdfPages):
~/miniconda3/lib/python3.7/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
39 renderer.start_filter()
40
---> 41 return draw(artist, renderer, *args, **kwargs)
42 finally:
43 if artist.get_agg_filter() is not None:
~/miniconda3/lib/python3.7/site-packages/matplotlib/figure.py in draw(self, renderer)
1862 self.patch.draw(renderer)
1863 mimage._draw_list_compositing_images(
-> 1864 renderer, self, artists, self.suppressComposite)
1865
1866 renderer.close_group('figure')
~/miniconda3/lib/python3.7/site-packages/matplotlib/image.py in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
129 if not_composite or not has_images:
130 for a in artists:
--> 131 a.draw(renderer)
132 else:
133 # Composite any adjacent images together
~/miniconda3/lib/python3.7/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
39 renderer.start_filter()
40
---> 41 return draw(artist, renderer, *args, **kwargs)
42 finally:
43 if artist.get_agg_filter() is not None:
~/miniconda3/lib/python3.7/site-packages/matplotlib/cbook/deprecation.py in wrapper(*inner_args, **inner_kwargs)
409 else deprecation_addendum,
410 **kwargs)
--> 411 return func(*inner_args, **inner_kwargs)
412
413 return wrapper
~/miniconda3/lib/python3.7/site-packages/matplotlib/axes/_base.py in draw(self, renderer, inframe)
2745 renderer.stop_rasterizing()
2746
-> 2747 mimage._draw_list_compositing_images(renderer, self, artists)
2748
2749 renderer.close_group('axes')
~/miniconda3/lib/python3.7/site-packages/matplotlib/image.py in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
129 if not_composite or not has_images:
130 for a in artists:
--> 131 a.draw(renderer)
132 else:
133 # Composite any adjacent images together
~/miniconda3/lib/python3.7/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
39 renderer.start_filter()
40
---> 41 return draw(artist, renderer, *args, **kwargs)
42 finally:
43 if artist.get_agg_filter() is not None:
~/miniconda3/lib/python3.7/site-packages/matplotlib/collections.py in draw(self, renderer)
410 self._linewidths, self._linestyles,
411 self._antialiaseds, self._urls,
--> 412 self._offset_position)
413
414 gc.restore()
~/miniconda3/lib/python3.7/site-packages/matplotlib/backends/backend_pdf.py in draw_path_collection(self, gc, master_transform, paths, all_transforms, offsets, offsetTrans, facecolors, edgecolors, linewidths, linestyles, antialiaseds, urls, offset_position)
1873 offsets, offsetTrans, facecolors, edgecolors,
1874 linewidths, linestyles, antialiaseds, urls,
-> 1875 offset_position)
1876
1877 padding = np.max(linewidths)
~/miniconda3/lib/python3.7/site-packages/matplotlib/backend_bases.py in draw_path_collection(self, gc, master_transform, paths, all_transforms, offsets, offsetTrans, facecolors, edgecolors, linewidths, linestyles, antialiaseds, urls, offset_position)
240 transform = transform.frozen()
241 transform.translate(xo, yo)
--> 242 self.draw_path(gc0, path, transform, rgbFace)
243
244 def draw_quad_mesh(self, gc, master_transform, meshWidth, meshHeight,
~/miniconda3/lib/python3.7/site-packages/matplotlib/backends/backend_pdf.py in draw_path(self, gc, path, transform, rgbFace)
1821 def draw_path(self, gc, path, transform, rgbFace=None):
1822 # docstring inherited
-> 1823 self.check_gc(gc, rgbFace)
1824 self.file.writePath(
1825 path, transform,
~/miniconda3/lib/python3.7/site-packages/matplotlib/backends/backend_pdf.py in check_gc(self, gc, fillcolor)
1771 delta = self.gc.delta(gc)
1772 if delta:
-> 1773 self.file.output(*delta)
1774
1775 # Restore gc to avoid unwanted side effects
~/miniconda3/lib/python3.7/site-packages/matplotlib/backends/backend_pdf.py in output(self, *data)
730
731 def output(self, *data):
--> 732 self.write(fill([pdfRepr(x) for x in data]))
733 self.write(b'\n')
734
~/miniconda3/lib/python3.7/site-packages/matplotlib/backends/backend_pdf.py in (.0)
730
731 def output(self, *data):
--> 732 self.write(fill([pdfRepr(x) for x in data]))
733 self.write(b'\n')
734
~/miniconda3/lib/python3.7/site-packages/matplotlib/backends/backend_pdf.py in pdfRepr(obj)
252 elif isinstance(obj, (float, np.floating)):
253 if not np.isfinite(obj):
--> 254 raise ValueError("Can only output finite numbers in PDF")
255 r = b"%.10f" % obj
256 return r.rstrip(b'0').rstrip(b'.')
ValueError: Can only output finite numbers in PDF

Add adapmap for dimensionality reduction

Hi @Xiaojieqiu ! I'm opening this issue to continue the conversation we started here.

I've made some small changes to utils_reduceDimension.py in order to add adapmap as a dimensionality reduction method (a non-linear alternative to PCA, essentially consisting of adaptive multiscale diffusion maps).

Layouts learned from this generalized diffusion structure are more informative than those learned on PCA embedding, and UMAP is particularly well-fit for embedding complex data due to a sequential approximation of the Laplace-Beltrami operator. Accordingly, dbMAP embeddings are able to visually depict cycling populations into cyclic trajectories and captures rich details from data structure into comprehensive embeddings.

The key idea is to explore visualizations derived from non-linear diffusion structure, rather than those derived from linear projections (PCA).

I'm new to contributing within GitHub, so any help to go through the pull request would be greatly appreciated. I had issues while following the contribution guidelines.

How much do landscape-derived attractor and transition pathway overlap with cell cluster and pseudo-time?

Thanks for developing this powerful software to analyze the dynamics of cells from single-cell data.

We recently used dynamo to analyze the CART cell clone dynamics from single-cell time course data. It helped us a lot to elucidate accelerated and diverted cell differentiation dynamics. However, there are several issues that we could like to discuss to enable more people to use this powerful tool for more convenient analysis.

It is cool to derive the pointed attractors and saddle points from RNA velocity data. From a theoretical point of view, pointed attractors should correspond to stable cell types and transition trajectories should correspond to single-cell transition manifold, which is usually inferred by pseudo-time.

How much do landscape-derived attractor and transition trajectories overlap with cell clusters and pseudo-time?

If we want to plot some specific gene’s RNA velocity, where can we find it ?

Usually, RNA velocity in 2D projection is decided by the top 30 ~ 50 genes.

However, since we are especially interested in certain predefined genes. If we want to plot these specific gene’s RNA velocity, where can we find it and how to implement it in the current APIs?